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PREFACE 


This report constitutes the second of two volumes which summarize 
the work accomplished under Contract NAS8-21282. It contains supporting 
experimental data, a theoretical analysis, and a listing of a digital computer 
program designed for predicting dynamic stability of propellant tanks under 
longitudinal excitation. The first part of the work, which deals with the 
dynamic state of the system prior to instability, is summarized in Final 
Report, Part I, entitled "Influence of a Rigid Top Mass on the Response of 
a Pressurized Cylinder Containing Liquid. " 

Both Part I and Part II of this Final Report are published on the 
same date. They present significant extensions and refinements of concepts 
originated under previous investigations conducted for NASA-MSFC. 

Results of this preliminary work are summarized in "Dynamic Stability and 
Parametric Resonance in Cylindrical Propellant Tanks, " by Daniel D. Kana, 
Wen-Hwa Chu, and Tom D. Dunham, Final Report, Contract No. NAS8- 
20329, SwRI Project No. 02-1876, January 17, 1968. 



ABSTRACT 


Dynamic instability and associated parametric resonance is a 
dominant form of response in a longitudinally excited cylindrical shell con- 
taining liquid. In order to assess the significance of such responses in a 
space vehicle propellant tank, the present paper is devoted to a theoretical 
and experimental study of their occurrence in a cylindrical shell system 
which includes the influences of axial preload, ullage pressure, partial 
liquid depth, and a finite top impedance. Donnell shell theory along with 
a modified Galerkin procedure is utilized to formulate equations which 
govern the stability of perturbations superimposed on an axisymmetric 
initial state of response. Stability boundaries are computed for a range 
of parameters affecting the region of principal parametric resonance and 
are compared with experimental results. It is found that liquid depth, top 
impedance, and ullage pressure have a strong influence on stability, while 
the effects of axial preload are relatively insignificant. 
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NOMENCLATURE 


a 

c 0 

c s 

E 

g 

H 

H s 

J z 
j t 


m 


N ixd* N led 

N 5xs» N 00s» 



radius of the shell 

speed of sound in. the liquid 

E/ p s , speed of stress waves in the shell 

modulus of elasticity 

standard acceleration of gravity 

h/ a, nondimensional liquid depth 

h s /a, nondimensional thickness of shell 

mass moment of inertia of top weight about z axis 

length of the shell 

one-half of the number of circumferential nodes; cos (m0) 

dynamic part of initial- state stress resultants [nondimen- 
sionalized by (1 - v^)/Eh g ] 

static part of initial-state stress resultants [nondimens ion - 
alized by (1 - v^)/Eh s l 


n 


axial wave number; sinnirx/i 


p r 

p o»Po 

R, B, X 
U, V, W 

X 0 

Z 0 

P 


nondimensional pressure loading on shell, P r /E 
axial preload, ullage pressure 

cylindrical coordinates (space-fixed) nondimen sionalized by 
radius a 

shell displacements u, v, w, nondimens ionalized by the 
radius a 

nondimensional amplitude of axial excitation (Xq = xq/ a) 
top acceleration impedance (force/ acceleration) 
density parameter p^a/p g h s 


v 



NOMENCLATURE (Cont'd) 


V 

Poissons ratio 

$ 

velocity potential, nondimensionalized by a)Qa /w r 

Pi 

mass density of liquid 

Ps 

mass density of the shell 

T 

nondimens ional time, T = co r t 

-1 

liquid parameter Cq 7 a^ 

co r 

response frequency 

0) 

excitation frequency 

<°k 

natural frequency of m-k'th mode 


designated frequency, nondimensionalized by a^/c| 

£2? 

X 

designated frequency, nondimensionalized by (1 - v^) X 

a 2 /c| 

n? 

~ i 

designated frequency, nondimensionalized by a ^/cq 

Superscripts 


n 

the amplitude of ( ) 

(•) 

(d/dT) ( ), T = wt 

( ) p 

related to initial state response 
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INTRODUCTION 


Dynamic instability and parametric resonance are known to occur in 
many engineering systems. Numerous classic examples have been studied 
in detail by Bolotin^. In more recent experimental investigations^, it was 
found that this type of behavior is dominant amidst a complex variety of 
responses which can be observed in a longitudinally excited model vehicle 
propellant tank which is not sufficiently reinforced with stiffeners. A theoret- 
ical and further experimental investigation^ was conducted for a longitu- 
dinally excited, liquid-filled cylindrical shell. It was found that the system 
initially tends to respond in a state comprised of linear axisymmetric modes. 
However, the resulting membrane stresses form a parametric load with 
respect to nonaxisymmetric perturbations superimposed on the initial state. 
Thus, for wide ranges of the excitation parameters, instability and subsequent 
parametric resonance results, and linear vibration theory is no longer ade- 
quate to predict the response of either liquid pressure or wall motion. 

In order to assess the significance of such instabilities in a propellant 
tank which forms a component in an overall space vehicle structure, the 
present paper is devoted to a study of their occurrence in a cylindrical shell 
system which includes the influences of axial preload, ullage pressure, 
partial liquid depth, and a finite top impedance. A diagram of the system is 
shown in Figure 1, which includes the appropriate parameters and boundary 
conditions for both the initial and perturbed states. The initial state 
represents linear forced axisymmetric motion, whose responses have already 
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Figure lb. Mechanism Of Dynamic Instability - 

Perturbed State 
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been determined, along with natural frequencies and modal functions for the 
system^. Stability of motion in the perturbed state is the subject of the 
present paper, although results from Reference 4 must be utilized in part 
of the analysis. Note that, theoretically, the perturbed state can be either 
axisymmetric (m = 0) or nonaxi symmetric (m > 0); however, for a single 
tank system, the nonaxisymmetric form of instability is dominant. 

DERIVATION OF STABILITY EQUATIONS 
The perturbed motion represented by Figure lb will be analyzed by 
means of Sander 's nonlinear shell equations 6 which are based on Donnell 
approximations. These equations contain nonlinear terms resulting from 
the rotation of shell elements as well as nonlinear strain-displacement 
relations. We will follow the philosophy of Bolotin* and assume that 
retaining only nonlinear terms which result from rotations is sufficient to 
determine dynamic stability. Compressible flow theory is used for the liquid. 

The motion is expanded into a series of the natural mode eigenvectors which 
were obtained from the solution of the free vibration problem^. A modified 
Galerkin procedure is then utilized to reduce the system to a linear second- 
order, time dependent set of coupled differential equations having periodic 
coefficients. The method is "modified" in the sense that the natural modal 
functions (finite series eigenvectors) are chosen as weighting functions, although 
they are not of closed form. An approximation of the perturbed motion will 
then be obtained by the use of only one eigenvector term of the series, so that 
the coupled set reduces to a single stability equation. 
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Thus, the governing shell equations are of the form 

F i = L U U + L 12 V + L 13 W - ^ 8 2 u/8t 2 _ o 

f 2 = l 21 u + l 22 v + l 23 w - Sja 2 y/ar 2 = o 

F 3 = L 3 iU + L 32 V + L 33 W - Q 2 9 2 W/ar 2 + e d L 33 W cos cot 


where 


_ a 2 , i - v a 2 

J 1 1 o n o 

ax 2 2 ae 2 


a 


J i3 " v ax 


T - li 4. 1 - v 9 2 

22 ^02 2 8 ^2 ' 


T2 


l 21 = 


"31 


"33 


ae* 


- V 


ax ’ 


l 32 - 


h|/_8_4 


+ 2 


a 4 


® 4 A 


+ i 


12 lax 4 ax 2 ae 2 ae 4 / 

+ ( Nx6s p) + fe ( N ® xs ?5c) + fe ( N ® 0S p) 


j.iL 

-v 2 ) 

T H s 

1 + V 

a 2 

2 

axae 

1 + V 

a 2 

2 

axae 

a 


ae 


a 


" ae 


I 


_ax ' 

^ xxs 

+ -( 
80 \ 

^ N ee s 


P = 0 
r 


t 9 / 3 \ 3 / ^ * 9 

L 33 -Qx[ K xxd axj + 8 0 V N 00d 90 


e d = 0 for free vibration 
e d = 1 for forced vibration 


(la) 

(lb) 


(lc) 


(2a) 


(2b) 


Pressure loading on the shell is given by: 
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(1 - v 2 ) - ?a 9$ 

T- p r = - 57 at R = 1 

s 

where the fluid velocity potential is governed by 

z 

V* $ - f2 2 9 2 $/9r 2 _ o 

Boundary conditions on the fluid and shell are: 

At X = 0, 

U = 0, W = 0, V = 0, 9 2 W/9X 2 = 0 

At X = l / a, 

W = 0, V = 0, 9 2 W/9X 2 = 0 

and 

f 4 = 9U/9X + z **n 2 a 2 u/ 9 r 2 = o 

where 


(3) 


(4) 


(5) 




Zp* 1 - 

2TT Ps a 2 h s 


for m = 0 


,** _ U - v 2 )Iz 

4p s h s a 5 


for m = 1 


Z** = oo for m > 2 


Solutions of the shell motion having a given circumferential displace- 
ment distribution will be sought as expansions of the m, k-th natural modes 


K K 

U(0,r,X) = cos m9 £ a k( T > U mk( X >= I a k U lk 

k= 1 k= 1 


K K 

V(e,T,X) = sin me £ a k (T)V mk (X)= £ a k U 2k 

k=l k = 1 


(6a) 


(6b) 
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K K 

W(0,t,X) = cos m6 £ a k (T)W mk (X) = X a k U 3k (6c) 

k = 1 k= 1 

where for convenience we have defined a general shell displacement 
vector 

U = U(U,V,W) = U(U 1 ,U 2 ,U 3 ) 

which is a function of both space and time and is associated with the fluid 
velocity potential $ and upper shell displacement . 

The potential $ satisfies Equation (4) which forms a constraint on 
the shell system. In order to interpret the fluid pressure loading as an 
apparent mass which is valid at the response frequency w r , we express the 
potential for forced motion as 

K 

$(0, T,o r , R, X) = cos m0 ]T a k (T)$ mk (u r> R, X) (7) 

k= 1 

Note that $ mk (w r , R, X) is the component of $ associated with a shell displace- 
ment component W mk , and both liquid and shell motion is anticipated to be 
nearly periodic with responses at frequency u r . At the shell wall, we use the 
notation 

For the special case of w r = w k , the system responds in the m, k-th 
natural mode, and the shell displacement modal functions form the vector 

Uk = U k (U lk , u 2k , U 3k ) 
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This vector is a function of space only and is associated with the fluid 
velocity potential X) and the top displacement Ug k . From the 

definition of natural frequencies, these modal functions satisfy 


3 

X ^ijUjk + X) + ^ k U| k = 0 

j = 1 

i = 1,2,3 


(8 a) 


au lk /ax - Z**n£u ik at X = i / a (8b) 

We now consider the forced motion. By means of a Galerkin 
procedure' , we form an expression for virtual work in the system 


I / Fi<U) • U ik , dS + E m / F 4 (U ik ) • U <k , d0 = 0 (9) 

i = 1 S C 


where 

e m = 1 for m = 0, 1 

= 0 for m > 2 
m — 

More specifically, we substitute Equations (3), (6) and (7) into Equations (1) 
and (5) and then by means of Equation (9) form an expression for virtual work 
between forces (expressed in terms of displacements) associated with the 
general forced motion and displacements associated with the m, k' -th natural 


mode. There results: 
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r » 3 3 

1 a k J I I Uik'< L ijUjk) 

= l L c i = l j = i 


k= 1 L s 


r ■ r 3 

- a k J X)U 3k . cos mQ dS - fi?a k J X U ik .U ik dS 

S S i = l 

f f / 9U lk 

+ € d a k cos wt ^ u 3k'( L 33 u 3k) dS + e m J (^a k U* k . 


~ / / 9U lk 

+ e d a k cos wt J U 3k .(L 33 U 3k ) dS + e m j (^a k U* k 


,) dS = 


+ Z"£i^ k U ik U; k , 


Upon use of Equations (8), reverting back to the more conventional displace- 
ment symbols in Equations (6) and carrying out the spatial integration, this 
can be written as 


X f ( K 2k'k - e m K 4k'kH^r a k + ^k a k> + K 3k'k(^r a k + ^ k a k )l 
k = 1 


+ X { a k^k^ M k'k^k^ + J k'k^ + a k^r^ M k'k^r^ + ■'■k'k^ 

k = 1 l 

+ a k N k i k cos ojtj- = 0 


^2k'k 


U mk' U mk dX 


i/a 

X 3k'k ~ J ^mk'^mkdX 
0 


K4k'k=Z**/ Uik'«| k de 
0 
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M k'kK> 


p 



/ 


0 


W mk - *mk(wr,X)dX 


(12d) 


u 2 ^ 

^k*k^ w k^ “ ^ ~~2 f (l^e) 

wf ./ 


i/a 

Vk = / W mk' W mk<iX (12£) 

0 


i/i 


N- 


k'k 


-/ 


w. 


mk 


TS.T* 

i v N. 


xxd 


d ^mk + ^xxd ®^mk 


8X‘ 


9X 


8X 


m ^00d W mk dX 


(12g) 


The coupled set of K equations (11) govern the perturbed motion, described 
in Figure lb, for a given value of m. We will limit further discussion to 
the case of modes having m > 2. For these modes, the dominant motion is 
radial for the set of natural modes at lower frequencies. Thus, the terms 
under the first summation can be neglected, and, in matrix notation, 
Equations (11) become 

6 2 [M] {a} + ?2 2 [K] |a} + [T] -[a} cos oat = 0 (13a) 

where the elements of the k-th row and k-th column of the corresponding 
matrices are 


¥k'k " M k'k^r) + Vk* 


^k’k " M k'k^ w k^ + X k'k 
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Tk*k “ ^k'k’ an< ^ ’ a r- ^ 


a l 

a 2 




a K 

Equation (13a) can further be written as 

Q Z 

{a} + ~ [Ml -1 [K] {aj- + [Ml -1 E X] {a} cos ut = 0 (13b) 

When the flow is incompressible, M^-i^- is independent of u and we have 
[Ml -1 [K] = [I] 

where [i] is the identity matrix. It must be emphasized that Equations (13) 
are not general differential equations in time but include the restriction of 
nearly periodic motion in the generalized apparent mass given by Equation (7). 


EVALUATION OF MATRIX ELEMENTS 

Modal Functions 

Elements of the matrices in Equations (13) will now be evaluated from 
Equations (12d-g) in terms of the X-dependent natural modal functions 
(eigenvectors) of the system. These functions, which are not of closed 
form, have previously been determined from an eigenvalue problem^ in 
terms of the following series forms for the shell displacements: 
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U mk< X > = \ B 2k x2 + B lk< X - 1 /») + B mOk 


+ 2 B mnk cos ^n X < 14a > 

n = 1 


V mk( X > = I C mnk sinX n X 
n~ 1 


(14b) 


W mk< X > =• I Amak sin \x X 
n= 1 


(14c) 


and, for the velocity potential. 


$ mkK* X ) = 2 A mnk*mnK’ X > 
n = 1 


where , X) is a component function which satisfies Equation (4) for 

vibration at frequency w r , and corresponds to the sin \ n X component function 
in W mk through the boundary condition which must be satisfied at the tank 
wall 4 . 

Mass Coefficients 

The mass coefficient I k i k will now be developed from Equation (12f). 
By means of Equation (14c), there results 


N N 


Z 2 A mnk A mn'k' s i n X n' X s * n X-n X ^ X 
n n = 1 n' = l 


N N ^ N 

7T 2 2 ^mnk^mn'k' ^n'n ~ TT 2 A mnk A mnk I 

<ia n = 1 n' = l n = 1 
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By substituting Equations (14c) and (15) into (12d), the 
liquid apparent mass coefficient corresponding to the response frequency 
co r becomes 

N N 

^k'k^r^ ^ ^ Vnk^mn'k^n'n^r^ ( 17 ) 

n = 1 n' = 1 

where 


01 


n'n^r ) 



$ mn (u r , X) sin \ n i.XdX 


Except for a normalizing constant, the latter expression has also been 
evaluated in previous work. That is, 


^n'nK> = a n' M : 


mn'n 


(18) 


where 


r i 

a n' = f sin 2 \ n iXdX = ^ 

0 

and M mn i n is given by Equation (18a) in Reference 4. Note, however, that the 
free index k used in the referenced expression is not the same k which is 
used to designate the natural mode herein, and we must also use w r = w. 

Finally, the apparent mass coefficient given by Equation (12e) 

is obtained simply by substituting a). r = in Equations (17) and (18). 
Parametric Coefficients 


Upon substitution of Equation (14c) into (12g), we obtain: 
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N N 

N k'k = X X A mnk A mn'k' J n'n (19a) 

n = 1 n' = 1 

where 


• / a / ** >{< 

Vn = / (xgNLd sin X n X - \ n cos \ n X 

0 

+ m 2 N| ed sin \ n X ^ s in \ n > X dX (19b) 

The dynamic stress resultant amplitudes N* xd and are produced by 

forced excitation in the axisymmetric initial state described in Figure la. 
These stress amplitudes can be expressed in terms of the amplitudes of the 
initial-state displacements by means of the usual stress-displacement 
equations 


NLd = ~^+ vW p 


(20a) 


Nee d = frP + v 


auP 

ax 


(20b) 


Thus, the parametric coefficients are partly determined by the initial state 
displacement amplitudes IjP and WP. 

The solution to the linear forced axisymmetric response of the initial 
state, in terms of displacements, has previously been given by Equation (25) 
in Reference 4. However, in the direct use of this equation, the appropriate 
elements of its matrices must include the substitution 

M** = Z** for m = 0, 1 


( 21 ) 
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since an arbitrary acceleration impedance is allowed in the present problem, 
rather than only a rigid mass. Further, to allow for comparison of numeri- 
cal and experimental data, it is convenient to express the initial-state 
displacements as ratios of the excitation amplitude Xq. Therefore, the 
dynamic displacement amplitudes are of the form 

5P = X 0 [i B|X* + Bf{x-i)+BP l0 

N 

+ X B mn" cos W' 
n" = 1 

N 

W P = X Q I A^" sin X n „X (22b) 

n" = 1 

whose coefficients are completely determined by solving for the case of 

x 0 = 1. 

The initial-state stresses can now be determined. Upon substituting 
Equations (22) into (20a), there results 

N -| 

= x 0 B? + B ? x - I (V'Bgin" - vA mn”) V' x (23a) 
L n" = 1 J 

and the derivative is 

3N* N *i 

-53T = X 0 [ B 2 - V' X 1 - vAgin") V' X J < 23b ) 

Upon substituting Equations (22) into (20b), there results: 
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Need = x o 


N 

v(B^ + B^X) + Y, (A^n" " v ^n" B mn n ) sin * 
n" = 1 


For convenience of computation, these stress resultants 
tives are expanded into complete Fourier series as follows: 

N 

= X 0 2 . N l n " sin K" X 

n" = 1 
N 

Ked = *0 2 N 2n „sinX n „X 

n" = 1 


3N: 


xxd 


N 


ax 


= Xq (B? + X N 3 n n cos X n "X 
\ 6 n" = 1 


where 


N ln" = B|X. ln „ + B?X 0n „ - + vA l 


mn 


N2n" ‘^'mn" v ^ B 2^1n" B 1^0n" “ V'®mn"^ 


N311" “ V'®mn" + v ^ , n"'^mn" 



X 


In" 


2a 

T~ 


SL /a 


/ 


i/a. 



s in X n tiX dX 


X sin \ n ,,XdX 


n-Xj (23c) 
and deriva- 


(24a) 


(24b) 


(24c) 


and 



The parametric coefficients can now be completely evaluated. 
Upon substitution of Equations (24) into Equations (19), there results 


N N 

^k'k " X 0 X X ^lnk'^nn'k' 


n = 1 n' = 1 


N 


X (Xn N ln" d n'Vn 


n" = l 


+ rc^N2 n nd n ii n , n - X n N2 n iie n n n i n ) - X n B2eQ n r n 


where 


(25) 


i/s 


: 0n 


i n = f cos \ n X sin \ n iXdX 


i/ a 

5 n i. n ' n = / cos X n „X cos X n X sin X n .XdX 


i/a 

d n"n'n = J sin X n " x sin X n ,x sin X n XdX 
0 

One -Term Approximation 

For the m-k'th natural mode in Equation (13a), set k' = k to obtain 


QyMajj + + Xj^Ta^ cos wt = 0 


(26) 


where from Equations (16-19) 


_ i 
M 


N N 


N 


- 2a X X XmnkVn'k^mn'n^r^ + 2a X \mk 
n - 1 n' = l n = 1 


TT t 


N N 


K 


* N 

& a 2 


2a 


X X + 2a ^ rnn ^ 


n = 1 n\ = 1 


n ■ 
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___ N N 

T = Z Z VnkVn'k 


N 

Z ^n^ln"^n"n'n. 


n = 1 n 1 = 1 


L n" = 1 


Z D 

+ m N2n n< ^n"n'n ” ^n^3n" e n"n'n) “ ^n®Z e On'n 

Equation (26) is a Mathieu equation whose stability properties are 
well known. To put it in a standard form® for determining the stability 
boundaries for principal parametric (1 / 2 -subharmonic) resonance, we set 
to = 2co r and obtain 


a^. + (a + 2c[Xq cos 2T)a^. = 0 

where 


(27) 


?2?K 

at ” 


— T 

q = 


2Sl*M 


The stability- boundaries can then be approximated by 


qXp = a - 1 for a > 1 


(28) 


qXg = 1 - a for a < 1 

In terms of input acceleration, which is convenient for experimental mea- 
surement, these become 

§x 



THEORETICAL AND EXPERIMENTAL RESULTS 
Experimental data for stability boundaries are obtained from the 
apparatus shown in Figure 2. All pertinent parameters, including the input 
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Figure 2. Experimental Apparatus 
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impedance Zq of the boundary condition at the tank top, could be measured. 
The use of acceleration impedance (force/acceleration) proved to be most 
convenient in this application. Variation of the impedance was achieved by 
using different rigid masses as well as the loading frame. The cylinder 
is made of 0. 005 -in. stainless steel, has a diameter of 10 in. , and is 
14.5 in. long (the same cylinder as that used in Reference 4). 

Theoretical and experimental stability boundaries are compared in 
Figure 3 for the k = 1, m = 10 mode. Theoretical results are obtained from 
Equation (29) with N = 5 terms. Excitation conditions at or above the 
boundaries result in a principal parametric resonance whose mode shape 
is dominantly the k = 1, m = 10 natural mode, and whose frequency of 
motion is 1 /2-subharmonic to the excitation. Experimental points were 
determined as the points of least acceleration where the parametric 
response would occur. It is apparent that significant deviation exists 
between theoretical and experimental results for the empty tank, and better 
agreement is achieved for greater liquid depths. After careful scrutiny, 
it was ascertained that the wider experimental stability boundaries are 
principally caused by imperfections in the cylinder. That is, split natural 
modes^ and spatially shifting modal patterns occurred so that one exact 
natural frequency did not exist. As a result, the experimental system shows 
a tendency to be more unstable than predicted by theory. This trend is 
apparent in all the data. It is possible that somewhat better agreement 
could be achieved by the use of some form of imperfection theory in the 
analysis. This possibility remains to be investigated. 



Input Acceleration , g x 





0 ' 0.994 0.996 0.998 1.000 1.002 1.0 

Frequency Parameter, w/2w,_ 10 
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Figure 3. Influence Of Liquid Depth On Stability 
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Figure 3. Influence Of Liquid Depth On Stability 
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It was desirable to determine the influence of the various system 
parameters on the stability boundaries for a given mode. This was done 
in terms of dimensional variables, in order to emphasize the complexity of 
this influence. For this purpose, it is necessary to understand the effects 
of the same parameters on the natural frequencies of the system. For 
convenience, some natural frequencies which were determined in the 
earlier work^ for several symmetric and one nonsymmetric mode are given 
as functions of liquid depth in Figure 4. 

It is recognized that, in general, a more unstable system will 
possess a stability boundary whose acceleration ordinate is at a lower value 
for a given value of the frequency parameter 2 <oj_jq. Therefore, in order 
to assess the effects of axial load, ullage pressure, liquid depth, and top 
impedance, a stability boundary acceleration g x j was determined at an 
excitation frequency value of = 0. 996 (Zojj.jq) f° r a range of each of 
these parameters. Theoretical and experimental results are compared in 
Figures 5 through 8. These results must be compared with those in Figure 4 
for proper interpretation. At a given liquid depth, increasing axial tension has 
only a small effect on natural frequencies and, likewise, only an insignificant 
effect on stability as shown in Figure 5. On the other hand, increasing ullage 
pressure significantly raises the natural frequencies of the nonsymmetric 
modes but leaves those of the lower symmetric modes essentially unchanged. 
Thus, as approaches a natural frequency for a symmetric mode, the 
parametric excitation of the initial state is amplified, and the system 
becomes more unstable. This is reflected by the dips in the curves in 
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Figure 4 . Natural Frequencies Of Partially Filled Tank 
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Figure 6. It is also interesting to note that, at certain frequencies, the 
system becomes completely stable where the parametric coefficient in 
Equation (27) becomes zero. 

Increasing liquid depth changes all natural frequencies, as shown 
in Figure 4, and has a profound influence on stability throughout the depth 
range, as shown in Figure 7. This results from the coincidence of w x ^ 
with natural frequencies of symmetric modes at certain points, as well as 
the provision of an increased distributed parametric loading on the tank 
wall. 

The influence of top impedance on stability is shown in Figure 8. 
Increasing this impedance lowers the frequencies of symmetric modes 
while leaving the nonsymmetric mode frequencies unaltered. Thus, strong 
interaction can again be seen to occur. The dip in the curve occurs at an 
impedance such that coincides with the natural frequency of the first 
symmetric mode. 

It is obvious that variation of the above parameters can cause either 
an increase or decrease of stability, depending on the range of analysis. 

Further, it must be recognized that many nonsymmetric modes are present 
in the frequency range indicated in Figure 4, and each mode can become 
unstable as the one which was studied. Therefore, a complex pattern of 
instability and parametric resonance occurs with many overlapping regions 
of instability. The overall trend of the data shows good qualitative agreement 
between theory and experiment, although significant quantitative discrepancies 
exist because of the reasons previously discussed. 






Input Acceleration, g X | 



CO 




Input Acceleration, g x , 



Top Impedance Z 0 , Ib/g 

Figure 8. Influence Of Top Impedance On Stability 
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INPUT DATA DESCRIPTION 


Card 

Fortran 

Variable t 



No. 

Symbol 

Name 

Units 

Definition 

1 

RHO 

386p 

lb/in. 3 

weight density of liquid 


RHOS 

386 Ps 

lb /in. 3 

weight density of the shell 


AO 

a 0 

in. 

inner radius of the tank 


SH 

h 

in. 

depth of liquid 


SHS 

h s 

in- 

thickness of shell 


SL 

l 

in. 

length of the shell 

2 

PO 

PO 

lb /in. 3 

ullage pressure 


ENU 

V 


Poissons ratio 


E 

E 

lb/in. 2 

modulus of elasticity 


BMSTR 

z** 


nondimensional top impedance 


CO 

c 0 

in. /sec 2 

speed of sound in the liquid 

3 

NJ 



no. of roots |i mj 


N 

N 


no. terms in series expressions 

4 

W 

co/2tr 

cps 

excitation frequency 


BM 

Z 0 

lb sec 2 /in. 

top impedance 


CPO 

p o 

lb 

applied force 


NOPT 



print option 


M 

m 


circumferential wave number 

5 

UMJ(I) 

^mj 


roots of (p- m j) =0, m * 0 

6 

W 

o)/2it 

cps 

excitation frequency 


BM 

Z 0 

lb sec^/ in. 

top impedance 


tNote that some variables in computer program are slightly different from 
those as defined in NOMENCLATURE on pp. v to vi of text portion of this 
report. 



Card 

Fortran 

Variable 



No. 

Symbol 

Name 

Units 

Definition 


CPO 

Po 

lb 

applied force 


NOPT 



print option 


M 

m 


circumferential wavenumber 

7 

UMJ(I) 

H-mj 

m = 0 

roots of (p m j ) = 0 , m = 0 


PROGRAM OUTPUT 

Printed Output 

1. Input data h, h g , f, ag, a, pg, Z**, 386p, 386p g , E, eg, v 


2 . 


3. 


4. 


In subroutine MITERS the mode no. , eigenvalue, no. of iterations, 
no. of times Aitken's delta process is used, the eigenvector, and 
check eigenvalue and eigenvector 


Mode no. , £2^, in rad/ sec, and co^ in cps 


co in cps, co^ in cps, K, T, M, a, q, Xq , g x> and co/2co^ 


PROGRAM NOTES 


Subprograms Used 

In addition to the main program the following subroutines are used: 

1. BES, computes the Bessel functions J n or X^-. 

2. MATINV, computes the inverse of a real matrix. 

3. MPRINT, prints matrix in matrix format. 

4. MITERS, computes the eigenvalues and eigenvectors of a real or 
complex matrix by the power method. 


5. 


NOPT is a print option that allows printing of intermediate results. 



The following subroutines are included as a part of subprogram MITERS 
SWEEPX 

NPNRMX 

DPMLTX 


SYMBOLIC LISTING 

Some of the program FORTRAN symbols which were not defined in 
the Input Data Description are: 


F ortran 
Symbol 

Variable t 
Name 


A 

a 


BMS 

M s = 2-irpg 

ah s i 

HS 

H s 


CS 

c s 


H 

H 


SLSTR 

f / a 


AONSQ 



XI 0 

X 10 


X20 

X 20 


ALN(I) 

X n 

n = 1,N 

ALNP(I) 

\ n ' 

n' = 1 , N 

XI (I) 

X ln' 

n' = 1, N 

X2(I) 

X 2n' 

n' = 1,N 

XOB(I) 

X 0n' 

n' = 1,N 

X1B(I) 

X ln' 

n 1 = 1, N 

ENB(I) 


n' = 1, N 


t Cor res ponding to Final Report Part I 



Fortran 

Symbol 

V ariable 
Name 



ETA(J) 

^mj 

j = l.NJ 


CJH(I) 

C jH 

j = l.NJ 


CNJ(I, J) 

Cn'j 

n’ = l.N; j = 

l.NJ 

BMJN(I, J) 

B . 
mjn 

j = l,NJ;n = 

1 » N 

BN1(I) 

^lmn 

n = 1, N 


BN2(I) 

N 2mn 

n = 1 ,rN 


BNO(I) 

^Omn 

n = 1 , N 


BOON(I) 

BQOn 

n = 1 , N 


BMMNN(I, J) 

■^mn'n 

n' = 1 , N; n = 

= 1,N 

R1M(I, J) 

^lmn'n 



S2M(I, J) 

^2mn'n 



U1M(I, J) 

U lmn’n 



U2M(1, J) 

^2mn'n 



T3M(I, J) 

^3mn'ti 



U3M(I, J) 

U3mn'n 



V2M(I, J) 

^2mn'n 



W1M(I, J) 

^lmn'n 



W2M(I, J) 

^2mn'n 



W3M(I, J) 

^3mn'n 



R2M(I, J) 

B-Zmti'n 



R3M(I, J) 

B-3mn'n 



S1M(1, J) 

^lmn'n 
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Fortran 

Symbol 

Variable 

Name 

S3M(I, J) 

^3mn'n 

T1M(I, J) 

^lmn'n 

T2M(I, J) 

^2mn'n 

V1M(I, J) 

^lmn'n 

X1M(I, 1) 

X, , 
In' 

U4M( 1 , 1) 

U 4 n 

V5M(1 , 1) 

V 5n 

R4M(1, 1) 

R 4n 

R5M(1 , 1) 

^-5n 

R6M( 1 , 1) 

®-6n 

01M(I, 1) 

°ln' 

PI M(I, 1) 

P ln' 

Y1M(1, 1) 

Y ln' 

Z1M(I, 1) 

Z In' 

Y2M(I, 1 ) 

Y 2n' 

Z2M(I, 1) 

Z 2n' 

Y3M(I, 1) 

Y 3n' 

Z3M(I, 1) 

Z 3n' 

P2M(I, 1 ) 

P 2n' 

Q2M(I, 1 ) 

Q2n' 

V4M(1 , 1) 

V 4n 

S5M(1 , 1) 

S 5 n 



Fortran 

Symbol 

Variable 

Name 

V3M(I, J) 

^3mn'n 

V5M(1, 1) 

V 5n 

S5M(1 , 1) 

S 5 n 

UVW(I, J) 

[ul 

RST(I, J) 

[R] 

UTR(I, J) 

[uJ-Hr] 

UWR(I, J) 

[[ U] - ?2 Z [ R ]] _1 

AI(I) 

I 0n' 

QBH(I) 

qB 

W 0n' 

FH(I) 

^rn' 

APH(I) 


BOPH 


B1PH 


B2PH 


CN1(I) 

N ln" 

CN2(I) 

N 2n" 

CN3(I) 

N 3n" 

DNNN(I, J,K) 

^n"n'n 

ENNN(I, J,K) 

e n"n'n 

E0NN(I, J) 

e 0n' n 

TPNN(I, J) 

T' . 
n n k 

BMBAR 

M 



Fortran 

Variable 

Symbol 

Name 

BKBAR 

K 

BTBAR 

T 

ABAR 

a 

QBAR 

TT 

XOAQ 

|(1- a)/q j 

XOSTR 

j(l - a)/ql 

WSTR 


WKBSQ 


WKSQ 

nZ 

“k 

AMNK(I, J) 

Vn'k 
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PROGRAM SLOSH ( INPUT, OUTPUT, TAPE60a INPUT ) 

determimation of natural frequencies and mode shapes 

of a partially filled shell 

PROJECT 02-1876 

SEPTEMBER 1967 

MODIFIED 4 SEPTEMBER 1968 

PROJECT 02-2332 

COMMON PI.H/COMl/A. AO,CSCOW/CO m2/UMJ( 5 T/C0M3/ AONSQ/COM4/ ALN ( 10 ) 1 

1 C0M5/aLNP( 10 ) 

DIMENSION XI < 10 ) , X2 1 10 > , X08( 10 ) » XlBMO ) » EN8(10),ETA< 5) 

1 ‘ ~ ' ‘ 

2 

3 

4 

5 

6 


C JH ( 5 ) , ON J ( 1 0 « . 5 ) , BM JN < 5 , 10 ) , BN0 { 10 ) , 8U ( 10 ) , BN2 (10) , 
BOON (1.0) ,BMMNN( 10,10 >,RlM( 10,10 ) ,R2M( 10,10 ),R3M< 10, 10), 
SIM < 10 * 10 ), S2 m< 10,10 >, S3MU0, 10 >, TIM <10,10 ),T2M< 10*10 >• 
t3m ( 10 , 10 ) , uIm ( 1 0 , 1 0 ) , u2m( 10 , 10 ) ,U3M(t0,10) # VlM(10,10), 
v2m<10,10), v3m< 10,10) ,WlM(10,10 ) , w2m(10,10),w3m(1O,10) , 

i Hi n < 7 7k ndw i 77 77 k ItUftMl 7 7 4 T i 4 n A A 4 Ilf OPIi , 77 7 A I 


UVW (33,33) 

dimens: 

l W5MC 
Z1MI 
Z3M I 

T 5 M t * , ->■>' ; , r,on4 > , *■ w , * 90 "i V ■* , * « / , I oru * * * I # u*n( * v , i , , r 4 - " » “ / I 

Q1M( 10, 1 ) # 02M ( 10 * 1) * P2M( 10 > 1 ) *Q2M ( 10 , 1) » 03M(10, 1 ) * P3M<10 » 1) * 
Q3M ( 10 , 1 ) , UTR < 33 , 3.6) , 8MMNNK ( 4,10,10) , BMMNN0 ( 4, 10 , 10 > 
DIMENSION IRON (34), IC0L<33), A 1 1 1 0 ) , QBW<10 ) , PH<33 ) , APH < 33 ) 
DIMENSION WCPS ( 4 ) , AMNK ( 33, 4 ) , CN1< 10) , CN2 ( 10 ) , CN3 <10 ) , 

1 DNNN{ 1 0 , 10 , 10 ) , ENNN ( 10 , 10 , 10 ) , E0NN(10 , 10 ) , TPNN< 10 , 10 ) 
DIMENSION GUESS! 33, 1 ) , VECTOR! 33, 4),EIGVAL( 4), 

1 NITER ( 4),US(33,8),HH(33,29),NAKSR< 4),EIGCPS( 5) 

PI A 3.14159265 
G a 386,0 
EPS = 1.0E-04 
SMLST s 1 , 0 E * 0 7 
NGUESS a 0 
NMODE = 1 
NITRSP a 100 
EPSP a 0 ,5E»°8 
AITKEN » ,9 
NTAPE a 1 

1000 READ 200, RHO,RHOS.a0,SH,SHS,SL 
200 FORMAT (7F10.0) 

IF ( EOF ,60)700,705 
STOP 

READ 205, P0,ENU,E.BMSTR,C0 
FORMAT <2F10,0,3E12.3) 

READ 215, N J, N 
FORMAT (315) 

N3 a 3 *n* 3 
NDIM a 33 
RHO « RMO/386, 

RHOS a RHOS/386. 

A s A0+Q . 5*SHS 
BMS a 2,0*PI*RH0S*A*SHS#SL 
HS 5 SHS/A 
CS a SQRTF<E/RH0S) 

H a SH/A 
SLSTR a SL/A 
AONSQ = 0,5*SLSTR 
X10 » 0,5*SLSTR 
X20 a SLSTR**2/3,0 
PRINT 505 

505 FORMAT <*1 PROJECT 02-2332*) 

PRINT 530, sh,shs,sl,ao, a,po,bmstr,rho,rhos,e,co,enu 


700 

705 

205 

215 


SL000100 
SL000200 
SL000300 
SL000400 
SL000500 
SL000600 
SL000700 
SL000800 
SL000900 
SL001000 
SL0 01100 
SL001200 
SL001300 
SL001400 
SL001500 
SL001600 
SL001700 
SL001800 
SL001900 
SL002000 
SL002100 
SL002200 
SL 00 23 00 
SL002400 
SL002500 
SL002600 
SL002700 
SL002800 
SL002900 
SL003000 
SL003100 
SL003200 
SL003300 
SL0034Q0 
SL003500 
SL003800 
SL003900 
SL004000 
SL004100 
SL004200 
SL004300 
SL004400 
SL004500 
SL004600 
SL0 0470 0 
SL004800 
SL004900 
SL005000 
SL005100 
SL005200 
SL005300 
SL005400 
SL005500 
SL005600 
SL005700 
SL005800 
SL005900 
SL 00 60 00 
SL0O6100 
SL006200 
SL0Q6300 
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530 FORMAT! *0 H * *.£10,3,* ( IN)*. 17X , *HS s *.£10,3, 

1* ( I M)* , 16X, *L * *,£10.3,* < !N)*//3X,*aO » *,E10.3,* <IN)*.18X, 
2 *A b *, Pi 0 .3, * (IN)*//* PO s *,E10.3,18H ( LB-SEC**?/ I N**4 ), 3x. 
36HM** = ,£10.3//* RHO * *,610,3.18H ( LB*SEC**2/ I N**4 ) , 2X » 

4 *RHOS b *,E10,3.18M <L8-SEC**2/IN**4)//4X. *E * *,E10.3, 

5 UN (LB/IN**2),11X,*C0 = *,£10.3,* U N/SEC >* , llX, *NU * *»E10,3> 
KOPT b o 

2000 READ 210, W , BM, cPO , NOPT , M 
210 FORMAT ( 3 F 10.0,215) 

READ 220, (UMJ<t),m,NJ) 

220 FORMAT (5F15.0) 

BM b BM/386. 

IF ( M ) 3 , 4 , 3 

4 BMSTR s ( BM* ( 1 , Q*£NU**2 ) ) / ( 2 , 0*PI*RHOS*A**3*KS ) 

JT PR s 1 

W « ,996*(2.*WCpS(JTER)> 
wRaD = 2 , *P I *W 
WK0 8 WcPS(JTER) 

WKSQ 8 ( 1 , -fcNU**2 ) * < 2 , *P I *WKB ) **2* ( A/CS )**2 
RO TO 2 
3 CONTINUE 

DO 7000 I TER*1 , 1 
IF (KOPT16000, 6005, 6000 
6000 w b .996*1 WCPSCITER) > 

WK8 8 WCPS(ITER) 

WKSQ s (1 ,-ENU**2)*(2.*PI*WKB)**2*< A/CS)**2 
6005 WRAD B 2 , *p I *w 

2 OMEGA * (WRAD*A)/CS 
WSO B omega ** 2 
WKBSO * (1.0*ENU*ENU)*MSQ 
CSCOW s (CS/CO )*OMEGA 
DO 10 I b 1 , N 
AIN(I) 8 ( J*PI*A)/Sl 
ALNP(I) » ( I *P I *A ) /SL 

Xl< I) s <2.0/(S| STR*ALNPC I) **2) )*(-!, 0*<-l)**t) 

x 2 ( n * < 4 , 0/ aln < 1 ) ** 2 ) * < «i )** 1 

XOB(I) B (2,0/(SLSTR*Ai.NP( n ) >* < 1 . 0- < -1 > ** I > 
xlp ( I ) 8 ( 2 , 0 * ( • 1 ) * * < I » 1 ) ) / A L NP ( I ) 

10 CONTINUE 
DO 5 1*1, N 
TMl e AL.NP < I ) +CSC0W 
TM2 * ALNPt I )-CSC0W 

ENB( I > * < ( <1 . O-COSF ( TMl*H ) ) /TMl ) * ( ( 1 , 0*COSF ( TM2*H ) ) /TM2 ) ) / ( 2 , 0* 
1 AONSO) 

5 CONTINUE 

DO 15 1*1, NJ 

E T A < I > * 50RTF ( aBSF (UMJC I )**2«cSc0w** 2) ) 
fF (UMJ( I )-CSCOw)20, 25,25 
25 TMl * EXPF(ETA( 1 )*M) 

CJH(I) * 0.5*(Tm1*(1.0/TM1)) 

60 TO 15 

20 CJH(I) 8 COSF(ETA( I)*H) 

15 CONTINUE 
DO 30 1*1, N 
TMl * ALNP< I >*H 
SN * SlNF(TMl) 

CN 8 COSE ( TMl ) 

DO 30 J*1,NJ 

IF (UMJ(J)-CSC0w>35,40,40 
40 TM2 a ETA<J)*H 
TM3 e EXPF(TM2) 

CSH b 0 ,5* (T*3*{ 1 . 0/TM3 ) ) 


SL006400 
SL 0 065 0 0 
SL006600 
SL006700 
SL006800 
SL006900 
SL007000 
SL007100 
SL007200 
SL 0 0 730 0 
S1007400 
SL007500 
SL007600 
SL007700 
St 00 79 0 0 
SL008000 
SL008100 
SL008200 
SL008300 
SU008400 
SL008500 
SLOO 86 00 
SL008700 
SL008800 

SL008900 

SL009000 

SLOO 9 IOO 

SL009200 

S1009300 

SL809400 

SL009700 

SL009600 

S1009900 

SL010000 

SlOlOlOO 

SL010200 

SL010300 

SL010400 

S1010500 

SL010600 

SL010700 

SL0108°0 

SL010900 

stonooo 

SlOlllOO 
SL011200 
SL011300 
SL011400 
SL 0 1150 0 
SL0H600 
SL011700 
SL01180Q 
SL011900 
SL012000 
SL012100 
SL012200 
SL012300 
SL012400 
SL012500 
SL012600 
SL012700 
SL012800 
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SNH a 0 i 5* < T^3“ ( 1 , Q/TM3 ) ) SL012900 
CNJ(I.J) a <M/ C aONSQ*(TM2**2+Tm 1**2) ))*<TMl*CN*CSH-TMl-TM2*sN*SNH) SLOl 3000 
GO TO 3H SLO 13100 


35 TH2 a ALNH< I)-Et.A(J) 

TK3 s AL^-Pf I)+ETA( J) 

CNJ< I, J) = (0 .5/A0NS0)*( ( <1,0-cOSF<Tm3*H> )/Tm3)*{ < 1 . 0-COSF < TM2* 

1 H))/TM2)> 

30 CONTI IVUE 

nO 45 1*1 , NJ 
DO 45 J* 1 , N 
SUMX a 0. 

SUMY * l,0£-50 
K s 0 

60 TRM a ( ( -1 )**K*pKN(K, J)*EMKJ(M,K, I > )/CJH< I ) 

SUMX = SUMX+TRM 


S1013200 
SLOT 3300 
SLO1340O 
SLOl 3500 
SL01 3600 
SlO 13700 
S1013800 
SL013900 
SLOl 4000 
Si. 01 41 0 0 
SLO14200 
SL01-43O0 


55 


50 

45 


IF ( ABSF < ( SUMX -SUM Y) /SUMY } -EPS >50*50,55 SL014400 

SUMY s SUMX SL014500 

K a K*1 SLOl 4600 

IF ( K * 1 0 0 ) 6 0 » 5 0 , 5 0 SL01470 0 

PMJN(1,J) a SUMy SL014800 

CONTINUE SLOI.49OO 

DO 65 1*1, N Si 0 1 5 0 0 0 

TMl * <1. 0-ENU*6NU)*WSQ SL015100 

TM2 a { a/SL)»0 ,5*BMSTR*RSQ SL015200 

TM3 b-(A/SL)+(-1)**I*BMSTR*WSQ SL015300 

TM4 a.{ a/SD*BMSTR*WSQ SL015400 

TM5 a 1 ,0-<A/SL)*Xl0 SL015500 

TM6 a 1,0*<M**2/(2,0*(1.0*ENU>*WSO) ) SL015600 

TM7 « 0 , 5* ( 1 , 0-pNU >*M*M*X20 SL015700 

RNO(I) a -((A/Si )**2*((2.0*TM1*X20-TM7)/(2.0*TM1))*(TM3/TM2)-TM5* SL015800 
1 TM6)/(1,0*( (2.0fTMl*X20-TM7)/(2. 0*TM1> ) * < TM4/TM2 > * < A/SL > **2- SLO 15900 

7 TM5*TM65 SL016000 

BNK!) s ♦< A/SL)*<1.0*BN0( m SL016100 

RN 2(1) a ( (TM4*BN0M)*TM3)*(A/SU**2)/TM2 SL01 6200 

CONTINUE SL016300 

DO 66 1*1, N SL016400 

SUMX = 0. SL016500 

SUMY a 1,06*50 SL016600 

K a 1 SLOl 67 00 

TRM a («1)**K*DKN(K, I ) *EMK j< M, K, 0 ) SL016800 

SUMX a SUMX*Trm SL016900 

IF ( ABSFf (SUMX-SUMY)/SUMY).EPS>67,67,68 SL017000 

SUMY a SUMX SLOl 71 00 

K a K*1 SL017200 

IF (K*l 0 0 ) 69# 67,67 SL017300 

BOON ( I ) a SUMY/COSF (CSC0W*H) SL017400 

CONTINUE SLOl 7500 

no 70 1*1, N SL017600 

DO 70 Js 1 , N SL017700 

SUMX a 0, SL017800 

SUMY a 1.0E-50 SLOI79OO 

IF ( M ) 86, 86 » 87 SLOlBOOO 

K a 1 SLOlglOO 

GO TO 85 SL01g200 

K* 0 SL 0 1830 0 

TRM ■ RMK(M,K)*DKN(K, J)*DKNB(K, I ) SL018400 

SUMX S SUMX+YRM SLO105OO 

IF (A6SF< (SUMX-SUMY) /SUMY) -EPS >75,75, 80 SL018600 

SUMY a SUMX SL01B700 

K a K*1 SL018800 

IF ( K-10 0 ) 85 1 85 , 75 SL018900 

75 SUMZ * 0. SLOI9OOO 


65 


69 


68 


67 

66 


86 

87 

85 


an 
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no 90 l ■» t # N J 

CALL besik, u*j<l>»o,xjo,t> 

TRH a XjO* 9 N*Ml. , J)*CNJ< I,L> 

SUMZ = SUM 2 +TRM 

90 CONTINUE 

If <M) 9 l. 9?,91 

9 ? CALL 8 FS(O,CSCOw,O*XJ 0 ,T) 

CALL BES(l,CSC 0 w,Q,XJl.T) 

TRM r -BO0N( J)*FNB< !)♦< ( ( 2. 0*6N8( I ) )/(CSCOw**2*COSF<CSCOw*H) ) )- 
1 ( (XJO*HKNB( 0. I ) )/<CSCOW*XJl) ) )*DKN(0. J) 

RNKNN(I,J) a ( (RhO*A)/(BWOS*SHS) )*(SUMY-SLMZ + TRm> 

GO TO 70 

9 1 TRM a 0 , 

BMKMN(I,J) a < (RH0*A)/(RH0S*SHS) >*(SUMY-SLMZ*TRM) 

IF ( ITER > 9 ^, 70,93 
93 CONTINUE 

IF (K 0 PT> 60 in, 6015,6010 
601 ft rHMNN 0 < ITER, I, J) a RMMNNU.J) 

GO TO 70 

6015 BMMNNK( ITER, l, J) a 0 MMNN( I , J) 

7 a CONTINUE 

IF ( N > 6025 , 6020,6025 
6025 CONTINUE 

IF <KOPT) 2000 , 6020, 2000 
6020 CONTINUE 

DO 95 I a 1 , N 
DO 95 Jat,N 
IF ( !-J) 100 , 105,100 
105 RlM< I , J) a 1 1 O + qMMNM I , J) 

GO Tn 95 

100 R 1 M( I , J) s BMMNNI I , J> 

95 CONTINUE 

DC 110 1 * 1 , N 
DO 110 J*l, N 
IF ( I»J> 1 15 , 120,115 
120 S 2 M( I, J) r 1 , 

GO TO nn 
115 S 2 m( I, J) a 0 , 

1 1 0 CONTINUE 

DO 125 lal.N 
DO 125 je 1 , N 

TNl a 1 , 0 +(HS*HS*(AIN( J)** 2 *M*M >** 2 >/ 12 , 

TN? a (P 0 /E)*( 0 , 5 *ALN{J)** 2 *M*M) 

IF ( NJ> 130 , 135,130 

135 X 3 N a ( A/SL>« ( , 5 *H** 2 -( ( 1 . -COS ( 2 . * ALN ( J)*H> )/( 4 .*ALN< J)** 2 > > ) 

TM 3 a X 3 N*H*M*(rhO*G*A/E) 

U 1 M ( I , J ) a m+(TM 2 *TM 3 -( ( CP 0 * ALN ( J ) **2 ) / < 2 , 0 *P I * A* A*E > > ) * ( < 1 . 0 - 
1 ENU**?)/HS> 

GO TO 125 

130 TM 3 « ALN(J>«ALNP(I> 

TM 4 a ALN( J)*ALNP( I ) 

X 3 N a ( A/SDM C ( 1 . -COS( TM 3 *M) >/Th 3 ** 2 ) -( ( 1 , -COS ( TM 4 *H > )/TM 4 ** 2 ) ) 
UIM(I.J) a (X 3 N*M*M>* ( ( 1 . 0 -ENU** 2 )/HS)*(RHG*G*A/E) 

125 CONTINUE 

no 140 I a 1 , N 
DO 140 Js 1 , N 
IF < l-J> 145 , 150,145 
150 u 2 m(I,J) s « 6 NU*ALN(J> 

T 3 M { I , J ) a 1,0 
U 3 M( I, JJ a M 

V 2 M { I , J ) a ALN ( J ) **2 + 0 . 5 * ( 1 , 0 « 6 NU ) *M*M 
W 1 M { I , J ) a M 


SL01 9100 
SL01 9200 
SL01 9300 
SL019400 
SLOl 9500 
SL019600 
SL019700 
SL019800 
SL01 9900 
SL020000 
SL 02 01 00 
SL020200 
SL020300 
SL020400 
SL020500 
SL02060 0 
SL020700 
SL020800 
SL 020 90 0 
SL021000 
SL0211O0 
Si 021200 
SLO2130O 
SL021400 
SL021500 
SL021600 
SL021700 
SL021800 
SL 021900 
SL022000 
SL022100 
SL022200 
GL022300 
SL 02240 0 
SL022500 
SL022600 
S1022700 
SL022800 
SL022900 
SL023000 
SL023100 
SL023200 
SL023300 
SL023400 
SL 02350 0 
SL023600 
SL023700 
SL 02380 0 
SL023900 
SL024000 
SL024100 
SL024200 
SL 02430 0 
SL024400 
SL024500 
SL024600 
SL024700 
SL024800 
SL024900 
SU025 000 
SL025100 
SL025200 
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W?M< !, J) x -0 . 5* ( 1 . 0*EN(J)*m*ALN<J) 

SL02530 0 


w3m<I,j) s 0 ,5*(1 . 0-FNU)*ALN< J)**2 + M+M 

SL025400 


<50 TO 140 

SL025500 

145 

U2KU,J) x Q, 

SL025600 


T.3m< I, j) x 0, 

Sl0?5700 


U3v< !, J) s 0, 

SL025800 


V2M( 1 , j) s 0 . 

SL 0 259 0 0 


WlM{ I , J) x 0 , 

SU026000 


U’2m { I , j) x 0 , 

SL026100 


W3M{ I , J) s 0 1 

SL026200 

14 0 

CONTINUE 

SL 02630 0 


PO 155 I si i N 

SL026400 


PC 155 Jsl.M 

SL 026500 


h2h(1,J) s 0, 

SL026600 


r3m( 1, j) s 0, 

SL 026 7 0 Q 


SIM ( I , J ) s 0. 

SL 02680 0 


S 3 m < ! , j ) s 0, 

SL026900 


TlM{J,J) b 0, 

SLO?7000 


T2M< J,J) s o, 

SL027100 

155 

CONTINUE 

SL027200 


PC 1 6 f' I x 1 » N 

Sl0?7300 


DO 160 Jx 1 , N 

SL027400 


IE U-J)l6l,l62 j l6l 

SL027500 

16? 

VlM< I , J) x -ALN( J)*ENU 

SL027600 


GO TP 160 

SL. 027700 

161 

v1m< I, J) x 0, 

SL027800 

160 

CONTINUE 

SL027900 


DO 165 I al > N 

SL028000 


X I'M (1,1 ) s X 2 M { f , 1 ) * X 3 M ( I » 1 ) a 0, 

Sl02q100 


u4m( 1,T) s u5M{1,I) * U6.M(1« I ) « 0. 

Sl0?8200 


V5M ( 1 , I ) x W 5 H < 1 1 I > = V 6H < 1 , I ) » W6M(1,I) * W4M(1,I) = 0. 

SL 02830 0 


B4m( 1 , i ) X S4M(1,n = T4M(l,n « 0, 

SL028400 


R5M(1, 1 ) X T5M(1, 1 J a 0. 

SL0285O0 


R6Y(1,I) a S6M(l,I> a T6M<1.1) a 0. 

SL028600 


Pi MU',!) x 0 ^ M ( 1,1) x 03 m< 1,1) a 0. 

SL 028700 


PlM (1,1) s P3M<!,1) b QlMn.l) a Q3 m<I,1) x 0, 

SL 028800 


Y1M< 1,1) s ENU*X08( I ) 

SL028900 


7lM{ 1,1) X 6NU*x1b( I ) 

SL029000 


Y 2 M ( 1 , 1 ) X M**2*< . 5*(1 . -ENU) )*xl< I > 

SL029100 


?2M<I,D s M**2*< .25*<1. -ENU) >*X2< I > 

SL029200 


y3m ( 1,1) X M*{ ,5+{l,+fcNj) )*X08( I) 

SL029300 


Z3m<!,i> x M*(,5#(l.+pNj))*Xl8<I) 

SL029400 


P2M< 1,1) X XI ( I ) 

S L 0 2 9 5 0 0 


<J2m ( f , 1 } X ,5*x?n> 

SL029600 


\/4 M ( 1 , I ) s -1 , 

SLO 297 OO 


S5M(l,l) s (*1)**I/(1.-ENU**2) 

SL029800 


no 165 Jxl , N 

SL029900 


IF (1* J)1 70»l75,l70 

SL 030 00 0 

175 

V3MU.J) x b ALN(J) *0.5*M*<1.0*6NU) 

SL030100 


GO TP 165 

SL030200 

170 

V3m ( I, j) x 0, 

SL030300 

165 

CONTINUE 

SL030400 


X5m a Z4M r 04M x P4M a Q4M x P5M a 0. 

SL030500 


X4M a -1, 

SL030600 


Y4M a SU/ A 

SL030700 


Y5M a 3 ,/BMSTR 

SL030800 


Z5M a Y4M*Y5M 

SL030900 


05M a 1 , / ( 1 , » E N l.J * * 2 ) 

SL031000 


05M a ( ,5MSI,/A)**2)/(1,-ENU**2) 

SL0311O0 


Y6M a ,5*(1.»ENU ))*<»(SL/A)*X10) 

SL 031200 


Z6M a -(l, r .5*<M**2*(l,-ENU> )*.5*X20) 

SL 0313 00 


X6m 0 , 5*m**2* ( 1 .-FNU) 

SL031400 
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06M a 1 , 

P6m a Xld-SL/A 
06M a . 5*X20 

IF <M)166,166, 16? 

1 67 X5M a 1 , 

Y5M a n, 

Z5M a. .b*<Sl/A>**2 
05 h a 0, 

P5M a 0. 

05M a 0, 

DO 168 1 a 1 , N 
V5M(1,I> a 
55H(1,!> a 0. 

168 CONTINUE 
166 HO 185 fal.N 

JK a UN 
JL a 1 ♦ ? * N 
IJ a 3*N+1 
IK a 3*N*2 
I L a 3 * N * 3 

UVW(I,IJ) = X 1 M ( 1,1) 
UVWU, IK) a Y.1M< 1,1) 

UV W ( I, ID a Z1M( 1,1) 
UVW < JK , I J ) 8 X2MU.1) 
UVWtJK, IK) a Y 2 w ( J , 1 ) 
UVW(JK.JL) a Z 2 m ( j , 1 ) 
UVW< JL, I J) a X3MI.1) 
UVW < JL > IK) a Y 3 m { 1,1) 
UVW < JL » ID S Z3m< 1,1) 
UVW ( I J, 1 ) a U 4 M ( 1 , I ) 

U V W ( I J , JK ) a V 4 M ( 1 , I ) 
UVW{ I J, JL > 8 W4M.fl, | ) 
UVW< IK, I ) 3 U 5 M ( 1 . i ) 

UVW { I-K , JK ) a V5M< 1,1) 
UVW < IK, JL ) a W5M ( 1, I ) 
UVW( IL, I > a U6M (1,1 ) 
UVW < IL , JK ) 8 V6 M < 1 , I ) 
UVWm .JL ) a W 6M < 1 , I > 
no 185 J a 1 , N 
IK a J*N 
IL a J*2«N 
UVW< I , J) a L1M( I, J) 
UVW( I, JK) a V 1 M ( I, J) 
UVW (I, ID a WlM<I,J) 
UVW ( JK, J) B U2M ( I , J ) 
UVW( JK, IK) a V2M ( I, J) 
UVW( JK, IL) = W?H( I, J) 
UVW( JL, J) * U3m < I , J) 
UVW(JL.IK) * V3M( I, J) 
UVW( JL, ID * WJM ( I, J) 
188 CONTINUE 

DO 190 I a 1 , N 
JK a I+N 
JL ■ !*2*N 
IJ a 3 * N* 1 
IK a 3*N*2 
IL a 3*N+3 

RST { I , I J ) a 01M( I ,1) 
RSTH.IK) a PlM(I.l) 
RST ( I, IL) = 01M( 1,1) 
RST( JK, I J) a 02 m( I , 1 ) 
RST( JK, IK) a P2M(I,1> 


SL '0315 0 0 

SL 031600 

Si 031700 

S L 0 3 1 8 0 0 

SL031900 

SL 0320 00 

S L 0 3 2 1 0 0 

SL032200 

SL032300 

SL032400 

SL032500 

SL 0 3260 0 

SL032700 

SL0.32800 

SL 0329 00 

SL 0330 00 

SL 033100 

SL03320 0 

SL033300 

SL033400 

SL033500 

S 1033600 

SL033700 

SL 033800 

SL033900 

SL034000 

SL034100 

SL034200 

SL 03430 0 

SL034400 

SL034500 

S L 0 3 4 6 0 0 

SL0347Q0 

SL034800 

SL034900 

SL035000 

SL 0351 00 

SL035200 

SL035300 

SL035400 

SL 03 55 00 

SL035600 

SL035700 

SL035800 

SL035900 

SL036000 

SL036100 

SL036200 

SL036300 

SL036400 

SL036500 

SL036600 

SL036700 

SL036800 

SL036900 

SL 037000 

SL037100 

SL037200 

SL037300 

SL 03 7400 

SL037500 

SLO 3760 0 
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PST< JK, It ) a GJ>m< 1,1) 

Si. 037 70 0 

BST(Jl.IJ) s 03* <1,1) 

SL 0 3 7 8 0 0 

RST< JL, JK) = P3m ( 1,1) 

SL03790D 

RSKJL.Il) = 03 m ( 1 ,1 ) 

SI. 038000 

RST ( J J, J ) a R 4 M ( 1,1 ) 

SL 0 * 8 1. 0 0 

RST< I J, JK ) a S4M(1, J ) 

SL03g2nO 

PST( I J , JL ) a T 4 M ( 1 , 1 ) 

SL 038300 

RST< IK, n a R5 m(1, 1 ) 

SL 0 3 8 4 0 0 

RST ( TK, JK) a S5 m < 1 , I ) 

S L 0 3 5 0 0 

RST (IK, JL ) a T5 m<1,I> 

SL 03 86 00 

RST( u , I ) a R6 m ( 1 , I ) 

SL 038 700 

RST I IL, JK) a S6 m < 1 , 1 ) 

S L 0 3 8 8 0 0 

RST ( I L * JL ) a T6 m<1,1) 

SL038900 

00 lpC Jri,N 

SL039000 

JK a J+N 

SL03q100 

IL a 

S L 0 3 9 2 0 0 

RST< 1 , J) a RlM ( f , J ) 

S L 0 3 9 3 0 0 

RST(l.lK) aSiHU.J) 

SL 0 394 0 0 

RST < I , IL) * TlM< I, J) 

SL 039500 

RST < JK , J) a W2M | I, J) 

SL039600 

RST< JK, JK) a S?M< 1, J) 

si 039700 

RST(JK,!L) a T2 m<1,J> 

SL039800 

RST< JL, J) a R3M ( 1 , J) 

SLO 399 OO 

RST( JL, IK) a S 3 m ( I, J) 

SL 04 00 00 

RSKJL.Il) a T3M ( I , J ) 

SL040100 

l9D CONTINUE 

SL040200 

1J a 3*N+1 

SL 040300 

IK a 3*N*2 

SL040400 

JL a 3 *n*3 

SL040500 

UVW< I J, I J) a X4 M 

SL040600 

UVW< I J. JK) a Y4M 

SL040700 

UVW ( I J, JL ) a 24 m 

SL040800 

U V W ( IK, J J ) a X5 m 

SL040900 

UVW ( IK, JK) a Y5 m 

SL041000 

UVW < IK, ID a 1\ j« 

SL041100 

UVW< JL , T J) a X6M 

SL041200 

uvw< IL, IK) a Y6M 

SL041300 

UVW( JL, ID « Z6* 

SL041400 

RST ( I J, J J ) a 04 m 

S1041500 

RST { I J, JK) a P4M 

SL041600 

RST < ! J, J.l ) a 04m 

SL041700 

RST ( IK, I J) a 05 m 

SL041800 

R S T ( IK, IK) a P5m 

SL 04 1900 

RST < I K , ID a Q5 m 

SL042000 

RST < IL , I J ) a 06 m 

SL042100 

RST ( IL, IK) a P6 m 

SL042200 

RST< IL, ID a 06 m 

SL042300 

IF ( NOPT ) 4 0 0 0 , 4n 05 , 40 0 0 

SL0424O0 

4000 PRINT 4010 

SL042500 

4010 PORMaT ( lHl , * < U > MATRIX *) 

SL042600 

CALL MPRINT ( U V w , N 3 , N 3 , n 0 I M ) 

SL 0427 00 

PRINT 4015 

SL042800 

4015 FORMAT (lHl,* < R ) MATRIX *) 

SL042900 

CALL MpRjNT < RST , N3 , N3 , NO I M ) 

SL043000 

4005 CONTINUE 

SLO<i3100 

IF ( M ) 198# *99, 19.8 

SL043200 

1 98 CONTINUE 

SL 04 33 00 

call matinv < UVW* I row , I COL * N3, NO I m » smlst > 

SL043400 

DO l9 ? I a 1 , N 3 

SL 04350 0 

00 1 9? Jb1,n3 

SL043600 

SUM a 0, 

SLO 43700 

DO l9l K a 1 , N 3 

SL043800 
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SUM s SUM * UVW( j ,K)*RST(K, J) 

191 CONTINUE 

UTRU,.i> = SUM 
UTRSV(I.j) s utr( i * j> 

19? cCNT I M ! E 

IF (NOPT)4920, 4025, 4020 
402ft PRINT 4030 

403ft FORMAT ( 1 H 1 , * ( U ) MATRIX INVERSE *) 

CALL MPRINT ( U V w « N 3 , N 3 , N D I M ) 

PRINT 4035 

4035 FORMAT ( 1 HI , « (u) INVERSE X (R) *) 

CALL MPRINT < JUT R , N3 , N 3 , N D I M ) 

4025 CONTINUE 

CALL MITERS (UTr,NTAPE,N3. GUESS, NGUESS.NMOCE, VECTOR, EIGVAL, 
1 N ITER , m I TRSP . EPSP , US , HH, M AXR , NC , A I TKFN, N AKSR , UTRSV ) 

DO 194 1=1,1 

XL DA = glGVALC I ). 

TM1 = (CS/A)/S0rT(1,-ENU**2) 

OMEGR s SORT (ABS( l./XLDA) ) *TM1 
CiMEGC S OMEGR/(?.*PI > 

IF ( ITER 1618,619,619 
619 E I GCPS ( I ) s CMEfiC 
618 PRINT 615, I, XI nA , OMEGR , OMEGC 
615 FORMAT ( 15,31:20 ,8) 

194 CONTINUE 

W » EIGCPSI ITER ) 

IF ( ITER >197,7000, 197 

197 WCPS(ITE») « <SoRT(APS(l./EIGVAL< ITER) ) >*TMl)/<2.*Pn 
CO 196 I «1 , N3 

AMNMI,ITER> b VECTOR ( I , ITER) 

196 CONTINUE 
7P0O CONTINUE 
KOPT b j. 

GO TO 3 
199 CONTINUE 

DO 195 Tsl,N3 
DO 195 Jsl , n3 

UWR ( I , J ) * UVW( j, J>-WKBS0*R5T( I, J) 

195 CONTINUE 

no 600 Isl,N 
J s N+T 
K s 2*N* I 

A I ( I ) s ( (-CSCOu*SINF ( AuNP( I )*H)*ALNP( I )*SINF(CSC0W*M> )/ 

1 (ALNPII >**2-CSCOW**2>)/AONSO 
QBH { I ) s (RHC*A*OMEGA*CO*Al<m/(RHOS*$HS*CS*COSF(CSCOW*H>) 
FH ( I ) b (1 .-6NU**2)*GBH< 1 > 

FH(U) s 0, 

FH(K> b 0, 

6&0 CONTINUE 

FH(3*N*1> a +1. 

FH(3*N*2) s 0 , 

FH ( 3*N*3 ) * 0 , 

IF (NOPTV4045. 4050, 4045 

4045 PRINT 4055 

4055 FORMAT { 1 hi , * ( u 5 - 0MEGA8SQ (R) *> 

CALL MPRINT ( UWR , N3, N3, ND I M ) 

PRINT 4 ft 6 0 

4060 FORMAT <*1 MATRIX <1)0 
PRINT 4065, (A I ( f ) , lsl,N) 

4065 FORMAT (5E20.8) 

PRINT 4070 

4070 format ( *0 matrix qhatb<i>*> 


SL043900 
SL044000 
SL0441O0 
SL044200 
Sl 044300 
SL044400 
SL 044500 
SL044600 
SL 044700 
SL044800 
SL 04490 0 
SL045000 
SL045100 
SL045200 
SL 0453 00 
SL045400 
SL045500 
SL 0 456 0 0 
SL045700 
SL045800 
SL 04590 0 
SL046000 
$L0461ft0 
SL046200 
SL046300 
SL046400 
SL 0 465 0 0 
SL046600 
SLP46700 
SL046800 
SL046900 
SL047000 
SL047100 
SL047200 
SL047300 
SL047400 
SL0A7500 
SL047600 
SL047700 
SL047800 
SL047900 
SL048000 
SL048100 
SL 04820 0 
SL048300 
SL048400 
SL048500 
SL048600 
SL048700 
SL 0 4880 0 
SL04890O 
SL049000 
SL049100 
SL049200 
SL049300 
SL049400 
SL 0 4 95 0 0 
SL049600 
SL 04970 0 
SL049800 
SL049900 
SL050000 
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PRINT 4065, (CRH( I ), 1=1, N) 

PRINT 4075 

4075 FORMAT ( * 0 MATRIX <F)*) 

PRINT 4065, (PHI I > , .1*1, N3) 

405n CONTINUE 

CALL PaTINV (UWR, IROW, 1C0L,N3,NDIM.SMLST) 

DO 605 I =1 , N3 
aPM< I > s 0 . 
nc 65 0 . j s 1 , N 3 

APH(I) a APH ( I ) ♦UWR ( 1 » J ) *FM( J ) 

61 0 CONTINUE 
6 O 5 CONTINUE 

IF (NOPTI4080, 4085, 4080 
4080 PRINT 4090 
4090 (-ORHAT (*0 APHAT(I) *> 

PRINT 4Q65, ( APw< I ), I *1 . N3 ) 

PRINT 4065, ( AlN< I ), Isl.N) 

PRINT 4065, (Xl( I). 1*1, N) 

PRINT 4065, ( X2 ( I ) , 1*1. N) 

PRINT 4065, (X0r( I), I sl.N) 

PRINT 4065, (XlR( !), Isl.N) 

PRINT 5040 

5040 FORMAT ( 1 M3. , * ( U ) - OMEGABSQ (R), INVERSE *) 

TALL MPRINT (UWR,N3,N3,NDIM) 

4085 CONTINUE 

pOpH s APH(3*N*1) 
pipH s APM(3*N+?) 

B2PH s APH( 3*N*3 ) 
no 2?b Isl.N 

J s M* I 

CM ( I ) s B2PR*XlPl( U*eiPH*XQ 8 < I > - A L N < I ) * APM «, ) + 6 MJ*APH < I ) 
CN? ( I ) s APH< n*ENtJ*(B2PH*xlBU )*81PE*X0B( I )-ALN< I )*aPH(v) ) 
C N 3 < I ) s ENU*At.N< I )*APH( I ) - aim ( I )**2*APH( J) 

225 CONTI-NUE 

no 230 Isl.N 
no 230 Jsl.N 
no 230 Ksi.N 
IF ( 1-J*K >805,800,805 

80n TERMl : n, 

TERMS = 0. 

TERM? s ( -1+ ( -1 )** ( J*K- I ) )/ ( aLN( J) + ALN( K) - AlN < I > ) 

TERM3 s f -1 ♦ f - 1 >**( )/(ALN(I )+ALN{ J)-aU(K) ) 

GO TO 830 

805 IF ( J*K-I >815,810,815 
8lo TERM? s 0, 

TERMl s (-1*(-1>**( J-J*K) )/tALNU)-ALN( J)*ALMK) ) . 

T£RM3 s (•!*( -!)**< I ♦ J ■- K ) ) / ( A L N ( I )*ALN(J)-ALN (K) ) 

TERMS s 1 1 , - ( - 1 ) * * ( J - I « K ) ) / ( A L N ( J ) - A L N ( I ) • A L N (X ) ) 

GO TO 830 

8l5 IE ( I * J»K )825, 820,825 
820 TERMS = 0. 

TERMl s (»1*<-1)**( I-J*K) )/(ALN( I )-ALN( J) + ALMK> ) 

TERM? s (-!■*( -1 )** ( J*K- I) )/( ALNf J) + ALN<K)-ALM I ) ) 

TERM5 s (1 .*(-1 )**( J-J-O )/( ALN( J) -ALN(1 l-ALMK) ) 

GO TO 830 

825 TERMl = (-1*(-1)**< 1-J*K> )/(ALN( 1 )-ALN( J)*ALMK> > 

TERM? s <-t+(-l>**< J+K-I ) )/< ALNC JUALNUKMALM I)) 

TERMS s (-1*(-1)**( I+J-K) )/CALN( I )+ALN( JMALMK) ) 

TERM5 s (3 ,-(-!)**( J- I-K) )/(ALN(J)-AL N( I )-Al.MX) ) 
p3o TERM4 = f *1»( -1 >**< I*j+<) )/( aLN( I )+ALN( J)*aLN (K) ) 

DNNN(I,J,K) s -f,25/A0NSQ)*(TERMl+TERM2+TEfiM3*TERM4> 
FNNN(I.J.K) s ( .25/A0NS5)*(TERm4*TpRm5-TER v 3-TFRm2) 


St 050100 
SI. 060200 
$1.060300 
SL05040 0 
SL060500 
SL060600 
SL 050 70 0 
SL050800 
SI. O 6 O 9 O 0 
SL051000 
SL061100 
SL051200 
SU051300 
St 061400 
SL061500 
SL051600 
SL 051700 

SL 061800 
SL051900 
SL052000 
SL0521 0 0 
SL052200 
SL062300 
SL062400 
SL052500 
SL 0526 0 0 
SL 0 62 7 0 0 
SL052800 
SLO 529 OO 
SL 05300 0 
S L 0 5 3 1 0 0 
SL053200 
SL0533O0 
SL053400 
SL053500 
SL053600 
SL053700 
SL 0538 0 0 
SL053900 
SL054000 
SL 0641 00 
SL 0542 0 0 
SL054300 
S1054400 
Si.05450 0 
SL 054600 
SL054700 
SL054800 
SL 06490 0 
SL 065000 
SL055100 

SL055200 
SL055300 
SL065400 
SL 0 555 0 0 
SL 055600 
SL 0 K 5 7 O 0 
SL055800 
SL055900 
SLP66000 
SL056100 
SL056200 
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230 CONTINUE 

DO 33 0 J s 1 , N 
00 310 K=1 , N 


SL 056300 
SL.056400 

SL 0565 0 0 


TERM1 s J*.K> )/<AUNC J)*ALN(K) ) SL056600 

IF ( J-K >840,835.840 SL0567OO 

835 TERM? r 0, SL056800 

GO TO 845 SL05690Q 

840 TERM2 s ( 1 , -<-!)** (J-K) ) / ( ALN< J) * ALN{ K ) ) SL057000 

845 E0NN<J,K) = ( ,!?5/A0NSQ> *< 2 , * ( TERM1 >TERM? ) ) SL057100 

310 CONTINUE SL057200 

M s 10 SL 057300 

DO 235 1=1, N SL057400 

DO 235 J = 1 » N SL057500 

SUM ■ 0, SL057600 

DO 240 K»1,N SL057700 

SUM = SUM* ( ALN ( J) #*2*CN1 < K)*DNNN < K, I , J) « ALN ( J) *CN3 (4 > *£NNN ( K , I « J ) SL0578QQ 
1 +M*M*nN2(K>*DNNN(K, I , J) ) SL057900 

240 cONTINUF SL05«000 


TPNN < I , J ) = CSUm-ALN< J)*B2PH*60NN( 1, J))*AMNK< J, JTER) 

235 CONTINUE 

IF (NOPT)4095, 5000,4095 

4095 PRINT 5005 

5005 FORMAT (1M1,12x,*N1(I)*,15X«*N2(I)*,15X»*N3(!)*) 

PRINT 5010, ( I , cNl ( I >,CN2(I> ,CN3( 1 > , I=1,N) 

5010 FORMAT ( I5.3E20 .8) 

PRINT 5035 

5035 FORMAT (*1 TP( I, J) O 

call mprint <tpnn,n,n,n> 

5000 SUmI = $UM2 = 0, 

DO 245 1=1, N 
SUM3 = SUM4 a 0 , 

DO 250 Jsl , N 

SUM3 = SUM3*AMNK< J, JTER)*8MMNN0( JTER, I, J) 

SUM4 = SUM4*AMNk( J, JTER)*BMMNNK< JTER. I. J> 

250 CONTINUE 

SUMl * SUMl*AMNK( I,JTER)**2*AMNK( I,JTER)*SUM3 
SUM2 « SUM2*AMNK(I,JTER)**2*AMNK(I,JTER)*SUM4 
245 CONTINUE 

RMBAR = WK8SQ*SuM1 
bkbar = WKS0*SUM2 

btbar = 0 , 

DO 255 1 = 1 , N 
DO 255 J«l, N 

BTRAR = RT8AR*TpNN( I , J)#AMNKM , JTER) 

255 CONTINUE 

ABAR = <4,*BKBAR)/RMBAR 

OBaR * (2,*bTbar>/bmbaR 

XOaQ = ABSF( <1 ,%ABAR>/QBAR> 

XOSTR = XOAQM <wRAD**2*A)/G) 

WSTR = W/(2,*WCPS( JTER)) 

PRINT 4040, 4, WKB, BKBAR, BTBAR, BMBAR-.ABAR, OBAR, XOAO, XOSTR, WSTR 
4040 FORMAT <1H0,2f10 .2.8E12.3) 

GO TO 1000 
END 

FUNCTION DKNIKD.ND) 

COMMON PI,H/COM4/ALN<10) 

IF (KD)10, 15,10 

15 DKN * (l.D-COSF( ALN(ND)*H) )/(H*ALN(ND)) 

GO TO 20 

10 DKN = ( ( 2 , 0*ALN(ND ) )/ <H* ( ALN(ND>**2* < (KD*PI l/H>**2) ) )*<!« 0” 

1 (-1 )**KD*COSF( ALN(ND)*H) ) 

20 RETURN 


SL058100 
SL058200 
SL058300 
SL058400 
SL058500 
SL058600 
SL058700 
SL058800 
SL058900 
SL059000 
SL059100 
SL059200 
SL059300 
SL05940Q 
SL059500 
SL059600 
SL059700 
SL059800 
SL059900 
SL060000 
SL060100 
SL060200 
SLOftOSOO 
SL060400 
SL06050 0 
SL060600 
SU060700 
SL060800 
SLO 6 O 9 OO 
SL061000 
SL061100 
SL061200 
SL06130Q 
SLO61400 
SL061600 
SL061700 
DK000100 
DK000200 
DK000300 
DK000400 
DK000500 
DK000600 
DK0Q070Q 
DK000800 
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EMD 

DK000900 


FUNCTION DKNB ( KnB , NPD > 

DK000100 


COMMON PI,H/COM3/AONSQ/COM5/ALNP(10) 

ok oo 02 no 


IF < KDB ) 10 , 15 *10 

DK000300 

is 

DKNB * (l,n'COSF(ALNP(NPD)*H>)/(AONSO*ALNP(NPD)) 

OK 0 0 0 4 0 0 


00 TO 20 

DK000500 

10 

DKNB = (ALNP(NPn)/AONSQ>*m,0-(-l)**KDB*CCSF(ALNP(NPD>*l-))/ 

DK000600 

1 

( AL NP ( NPD)**2»<.(KDB*P!) /H ) **2 ) ) 

OK 0 0 0 7 0 0 

20 

RETURN 

OK 0 0 0 80 0 


END 

OKOOO 9 OO 


FUNCTION RMK (Mr,KR) 

RM000100 


COMMON P1,H/COM1/a,aO,CSCQW 

RM000200 


DIMENSION T C 1 0 0 n > 

RM000300 


X I K * SQRTFt ABSFI ( (KR*PI*AO>/(N*A) )**2-CSC0W*#2) ) 

RMQ00400 


IF <KR-(CScOw*H)/P! >10, 10,15 

RM000500 

In 

CALL BPS ( MR, X I K , 0 , R JM, T ) 

RM000600 


CALL BES (MR+1,xIK,0»«JM1,T> 

RM000700 


RMK • RJM/(*MR*RJM-XIK*RJM1> 

RM000800 


BO TO ?Q 

RM 0 0 09 0 0 

15 

CALL BBS (MR.XlK.l.RIM.T) 

RM001000 


CALL E»FS (MR*1,xIK,1»R!m1,t> 

RM001100 


RMK B RIM/(t.MR*RIMfXlK*RIMl> 

RM001200 

20 

RETURN 

RM001300 


END 

RM001400 


FUNCTION EMKJ (mE,KE,JE) 

EM000100 


COMMON PI* H/C0M1/ A » A O » CSC0w/COM2/UM J ( ■5) 

EM000200 


DIMENSION T(10Pn> 

BM 0 0 0 3 0 0 


CALL BES (ME.UMjt JE>,O.EJM,T> 

6M000400 


ALFSO b 0,6*(l,n.(ME**2/UMJ(JE)**2))*f:jM**2 

EM000500 


X I K S SQRTF< ABSFI ( (KE*PI*aO>/(H*A) >**2-CSC0W**2 > ) 

EM000600 


IF < KE- < CSCOR*H ) /p 1 ) 10 » 15* 15 

EM000700 

In 

EMKJ b ,FJM/< (XIK**2-UMJ( JE)**2)*ALFS0> 

FMOOOaOO 


00 TO 20 

FM000900 

15 

EMKJ » +EJM/( (XJK**2+UMJ( JE)**2)*ALFSQ) 

EM001000 

20 

return 

EMOOllOO 


end 

EM001200 


subroutine be s< no, x,kode, result, T) 

BES 1 

03 UcSD BES 8ESS5L FUNCTION 

F-63 

c 

C3 UCSD BES 

E 63 


DIMENSION T < 1 0 0 o » 

BES 2 

107 

FORMAT (55H1 NEGATIVE ORDER NOT ACCEPTED IN BESSEL FUNCTION ROUTINE)0ES 3 


KLAMbI 

BES 4 


KQsNO*l 

BES 5 


IF ( X ) 6,1,6 

BES 6 

1 

I F ( NO ) 4,2,3 

BES 7 

2 

T ( KO > =1 , 0 

BES 8 


RESULTsl.O 

BES 9 


RETURN 

BES 10 

4 

IF (KO) 5,10,3 

8ES 11 

3 

RESULTsO 

BES 12 


RETURN 

BES 13 

10 

RESULTS?. 9R9999999E200 

BES 14 


RETURN 

BES 15 

5 

PRINT 107 

BES 16 


STOPl 

BES 17 

6 

IF(NO) 5,7,7 

BES 18 

7 

! F ( KODE ) 8,9,8 

BES 19 

8 

KLAM*KLAM*1 

BES 20 

9 

j0b2*1FIX(X) 

BES 21 


MOsNO 

BBS ?2 


IF ( MO*' JO > 11,12,12 

BES 23 

11 

M0 = JO 

RES 24 
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12 M0sM0*tl 

T(MO)=0. 

I. URsMO-1 
T ( LUB > *1 . QS-25Q 
GO TO (23, 51), K|. AM 

23 F =2*1 UR 
MO=MO-3 
I 2s MO 

24 F s F » 2 , 

T( 12*1) =F/X*T( 12*2 )*T< 12*3) 

IF < 12)25,26.25 

25 12=12-1 
GO TO 24 

26 SUMsT(l) 

1)040 Jr 3, MO, 2 

40 SUMsSUM+2 . *T ( J ) 

F«1 ./SUM 
DO 50 Js 1 , KO 

50 T(J)=T(J)*F 
RESUlTrT(KO) 

RETURN 

51 F»2*tUR.2 
M0aM0-3 
r2 = M0 

511 T ( I2*1)=F/X*T< 12*2) *T <12*3 ) 

I F ( 12)52,53,52 

52 1 2= 1 2-1 
F=F»2, 

GO TO 511 

53 SUMpT(I) . 

DO 70 J*2, MO 
70 SUM«SUM*2,*T(J) 

f=i , /Sum* E xPF ( x ) 

DO 80 J* 1 » KO 
80 T ( J ) «T ( J ) *F 
R6SULT=T(K0) 

RETURN 
END 

SURROUTINE MaTTnV ( A , I ROW i ICOL , N , NDIM , SMLST ) 
DIMENSION A ( 1 > , I ROW < 1 ) , I COL < 1 > 

709-16065 

709-16065 SUBROUTINE MATINV - MATRIX INVERSION ROUTINE 

a = array name of matrix 
I ROW s DIMENSIONED at n*i or greater 
ICOL = DIMENSIONED AT N OR GREATER 
N s NUMBER OF equations 

NDIM s VALUE OF I IN DIMENSION A(I,J) , l AND J MAY DIFFER 
SMLST 8 SMALLEST LEADING ELEMENT ALLOWED BEFORE CALLING THE 
SYSTEM SINGULAR » USUALLY s 1.0 6-04 OR 1.0 E-05 

NPl « N * 1 
DO 100 I b 1, m 
ICOL ( I ) = I 
100 IROW ( I > = I 

DO 240 ITER = 1, N 
MAXR « ITER 
MAXC » 1 

TEMP = aBSF ( A < MAXR ) ) 

LIMITC ■ NPl - I TRR 
DO 120 I = ITER, N 
DO 120 J = 1, LIMITC 


RES 

25 

RES 

26 

BfcS 

2 7 

8ES 

28 

RES 

29 

RES 

30 

BES 

31 

BES 

32 

BES 

33 

8ES 

34 

RES 

35 

BES 

36 

BES 

37 

RES 

38 

RES 

39 

8ES 

40 

BES 

41 

BES 

42 

BES 

43 

RES 

44 

8ES 

45 

BES 

46 

RES 

4 7 

BES 

48 

RES 

49 

BES 

50 

BES 

51 

BES 

52 

BES 

53 

BES 

54 

BES 

55 

0E 5 

56 

BES 

57 

BES 

58 

BES 

59 

9ES 

60 

BES 

61 

BES 

62 

MA 

100 

HA 

20 0 

HA 

300 

HA 

400 

MA 

500 

ma 

600 

MA 

70 0 

HA 

80 0 

ma 

9 0 0 

ha 

1000 

MA 

1100 

ma 

1200 

HA 

1300 

MA 

140 0 

MA 

1500 

HA 

1600 

MA 

1700 

MA 

1800 

ma 

1900 

MA 

2000 

MA 

2100 

MA 

2200 

MA 

2300 

MA 

2400 
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1 J * ( J - 1 ) • MD I M + I 

IF ( tfmp - ( AR$F ( A ( 1 J ) ) ) ) 110, 120, 120 

C 

llo MAXR s I 
MAXC = I 

TEMP s AHSF ( A ( 1 J ) ) 

120 continue 

IF ( TfmP - SML.ST ) 130, 130, 140 

c 

130 IROW ( NPl ) a JTER 
PRINT 1, ITgR , SMLST 

return 

C 

140 IF ( MAXR - ITER ) 150, 170, 150 

C 

150 no 160 J s 1, i\j 

MAXRj a ( J v 1 ) * NO 1 M ♦ MAXR 
ITJ s < J - 1 > * NDIM ♦ ITER 
TEMP s a < M A X R J ) 

A < HAXRJ ) s A < IT.J ) 

1 60 A t IT.. I ) a TEMP 

ITEMP s IRnw ( MAXR ) 

I Row < MAXR ) 8 I ROW ( I TER > 

TROW ( ITER ) s ITEMP 
170 IF ( MAXC - 1 ) lflO, 200, ,.l.gO 

c 

180 DO 190 I si, v 

1MAXC 8 < M A X C - 1 ) * AID I M * J 
TEMP s a ( 1 ) 

A ( I ) s A ( I MAXC ) 

190 A ( IMAXC ) a TRMP 

ITEMP = } COL ( M A X C ) 

I COL ( MAXC > a I COL ( 1 ) 

TOOL ( 1 > t ITfMP 

200 TEMP i A ( ITER ) 

ITEMP = ICOL < 1 ) 

DO 210 J s 2, m 
ITJM1 s ( J * 2 ) * NDIM * ITER 

I T J * < J - 1 > * NDIM ♦ ITER 

A < ITJMI ) 8 a < I T J ) / TEMP 
210 ICOL ( J - 1 >s ICOL ( J > 

JTN s ( N - 1 ) * NDIM * ITER 

A ( ITN ) s 1 . 0 / TEMP 

ICOL < N > a ITpMP 
DO 240 I a 1, M 
IF ( I « ITER ) 220, 240, 220 

C 

220 TEMP a a ( ! ) 

DO 230 j 8 2, m 
IJM 1 a ( J - 2 ) * NDIM ♦ I 
IJ a ( J - 1 > * NDIM * I 
ITjMl s ( J * 2 ) * NDIM ♦ ITER 
A ( IJMi ) s A ( u ) - A ( I T.JM.1 1 * TEMP 
230 CONTINUE 

IN s { N - 1 ) * NDIM + I 

ITN a ( N - 1 ) * NDIM ♦ ITER 

A ! IN ) a - < TEMP * A ( ITN ) ) 

240 CONTINUE 

DO 290 I = 1, n 
no ?50 J s I, n 

IF ( I ROW ( v ) - I ) 250, 260, 250 


MA 2500 
MA 2600 
MA 2700 
MA 2800 
MA 2R0 0 
MA 3000 
MA 3100 
MA 3200 
MA 3300 
MA 3400 
MA 3500 
MA 3600 
MA 3700 
MA 3800 
MA 3900 
Ma 4000 
MA 4100 
MA 4200 
M A 430 0 
MA 4400 
ma 4500 
MA 4600 
MA 4700 
MA 4800 
MA 4900 
MA 5000 
MA 5100 
MA 5200 
MA 5300 
MA 5400 
MA 5500 
MA 5600 
MA 5700 
MA 5800 
MA 5900 
MA 6000 
ma 6100 
MA 6200 
MA 6300 
MA 6400 
MA 6500 
MA 6600 
MA 6700 
MA 6800 
MA 6900 
M A 7 0 0 0 
MA 7100 
MA 7200 
MA 7300 
MA 7400 
MA 7500 
MA 7600 
MA 7700 
ma 7800 
ma 7900 
MA 80 0 0 

ma 8ino 
Ma 8200 
ma 83 0 0 
Ma 8*00 
Ma 8500 
MA 8600 
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25n 

CONTINUE 


HA 

870 0 

?6d 

IE ( 1 * J ) 270, 290, 270 


HA 

88 0 0 




HA 

8900 

270 

no 280 LSI, n 


HA 

9000 


1.1 S ( I - 1 ) * NDIM + L 


HA 

9100 


1 J 3 ( J - 1 > # ndim * L 


HA 

920 0 


TEMP = A < LI ) 


HA 

9300 


A ( LI ) s A { LJ ) 


HA 

9400 

28f> 

A ( LJ > = TEMP 


HA 

9500 


I ROW ( J ) s I Row < I ) 


MA 

9600 

291 

CONTINUE 


HA 

9700 


no 340 I = 1, M 


HA 

9800 


DO 300 J s I, m 


MA 

9900 


If « IrnL ( J ) * J ) 300, 310, 300 


ha 

10000 




ha 

10100 

300 

CONTINUE 


MA 

1020 0 

31 0 

If < I * J ) 3?0, 340, 320 


ma 

1 030 0 




MA 

10400 

320 

DO 330 L = 1, N 


HA 

10500 


IL s < L - 1 ) * NDIM + I 


MA 

10600 


JL s ( L - 1 ) * NDIM •* J 


ma 

10700 


TEMP 3 a ( IL ) 


HA 

10800 


A < IL ) s a ( JL ) 


Ha 

IO 9 OO 

330 

A ( Jl ) s TEMP 


HA 

110 0 0 


1 COL < J > * I COL ( I > 


MA 

11100 

340 

CONTINUE 


ma 

11200 


IROW ( NPl ) s 0 


Ha 

11300 


RETURN 


MA 

11400 




MA 

11500 

1 

FORMAT (7H0ON TWEI3.63HTH ITERATION ALL 

THE REMAINING 

TERNS WERE LMA 

11600 

*ESS THAN! OR EOUaL TO Ell. 4,18* IN ABSOLUTE VALUE) 

MA 

11700 


END 


MA 

11000 


SUBROUTINE MPRInT ( A , M , N , MO ) 


MP 

10 0 


MATRIX PRINT SUBROUTINE 


MR 

200 


THE CALL FOR TH t S SUBROUTINE IS AS FOLLOWS, 

MR 

300 


CALL MPRINT (A,M,N»MD) 


MR 

40 0 


WHERE A IS THe MATRIX TO BE PRINTED 


MR 

50 0 


M IS THp number of rows 


MR 

60 0 


N IS THf NUMBER OF COLUMNS 


MP 

7 OO 


MD IS DIMENSIONED NO, OF ROWS OF 

MATRIX A 

NR 

800 


DIMENSION A ( 1 ) , JT ( 6 ) , C ( 6 ) 


MR 

900 


EQUIVALENCE < JT , C ) 


MP 

10 0 0 


Ml 3 M 


MR 

1100 


n2 s 6 


MP 

1200 


n3 b 6 


HR 

1300 


N4 b 1 


MR 

1400 

100 

IF < N3 - Nl ) 120, 120, 110 


MP 

1500 




MR 

1600 

110 

N2 s Nl - N3 ♦ 6 


MP 

1700 


N3 b Nl 


MP 

leoo 

I2n 

K s 0 


HP 

1900 


DO 130 I s N4, N3 


MR 

2000 


K s K ♦ 1 


HP 

2100 

130 

JT I K ) s ! 


MP 

2200 


PRINT 1, ( -T ( I ) , I b 1 , N2 ) 


MR 

2300 


no 150 1 s 1, m 


MR 

2400 


K 3 0 


MR 

2500 


L 3 MD * ( N4 - 1 ) * I 


MP 

2600 


no 140 J 3 N4, N3 


MP 

2700 


K s K + 1 


MP 

2000 


C ( K ) b A ( L ) 


MR 

2900 

140 

L s L * MD 


MR 

3000 
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150 PRINT ?, I , ( C < K ) , K * 1 , N2 ) 

IT ( n3 - N 1 ) 160, 170, 178 

C 

160 N3 s N3 * 6 
n 4 s M » fi 
fi 0 TO 100 
C 

170 RETURN 
C 

1 FORMAT <1H , 4X, 6( 6X » 7HCQLUMN 1)4 ) / ) 

? FORMAT (1H 114, X, <6E 17,8) ) 

FNO 

SURRCUT I NE NPNRmX (A, 8, N, FI, INDEX, MD, NX ) 

C CALL NPNRMX, (A, B, N, Ft, INDEX, MD, NX, NP ) 

C AsVECT HR TO 0E FORMALIZED B»NORM AL I ZED VECTOR (May. a) 

C NsSIZE FUNORMaLIZINO NUMBER 

C INDEX** ON ENTRY, NORMALIZE ON NUMBER WHOSE INDEX IS INDEX 

C *0 ON ENTRY, NORMALIZE ON LARGEST S,H. AND SET 

C INDEX* TO ITS INDEX, 

C »- ON ENTRY, NORMALIZE ON Ft, 

C MDsSINGLE PRECISION DIMENSIONED NUMBER OF ROWS OF a AND B 

C NXel, VECTOR REaL 

C =2, VECTOR COMPLEX 

C NP si, SINGLE PRECISION 

DIMENSION A ( 1 ) , 8(D, Ft (1 ) , D(i), CM) 

N1 b 2 

k2*N 

N4bMD 

IF < INDEX) 32, 7,38 
7 GOTO <11, 8), NX 
A FLs < A(1)**2 *A(n4+1)**2> 

INDEXbI 
DO 10 K a N 1 , N 2 
I «K + N4 

D* •( A ( K ) **2 + A ( n**2) 

IF (FL-D > 9,9,10 

9 FI =D 
INDEXbK 

10 CONTINUE 
6 FLbA( INDEX) 

GOTO 18 

11 FL«ABSF(A(1).) 
iNDEXsi 

DO 13 KbN1,N2 
D« ABSF ( A < K) > 

IF (FL-D) 12,12,13 

12 FL»D 
INDEXbK 

13 CONTINUE 

14 FLsA( INDEX) 

16 DO 17 1*1, N 

17 B<I)bA(I)/FL 
GOTO 30 

lg I»INDEX*MD 
FL(2 )=a( I) 

1 9 n«FL(l)**2+FL(2)**2 

no 20 1 a 1 , N 

K«!*MD 

C» A ( I ) *FL ( 2 ) - a ( K ) *FL ( 1 ) 

r( n*(A( i>*fl(D*a(k)*fl(2))/d 

20 B ( K ) auC/D 
30 RETURN 


MP 3100 

MP 3200 

MP 3300 

MP 3400 

MP 3500 

MP 3600 

MP 3700 

MP 3800 

MP 3$00 

MP 4000 

MP 4100 

MP 4200 

MTRS0148 

MTRS0135 

MTRS0136 

MTRS0137 

MTRS0138 

MTRS0139 

MTRS0140 

MTRS0141 

MTRS0142 

MTRS0143 

MTRS0144 

MTRS0145 

MTRS0150 

MTRS0152 

MTRS0153 

MTRS0154 

MTRS0155 

MTRS0156 

MTRS0157 

MTRS0158 

MTRS0259 

MTRS0160 

MTRS0161 

MTRS0162 

MTRS0163 

MTRS0164 

MTRS0165 

MTRS0166 

MTRS0167 

MTRS0169 

MTRS0170 

MTRS0171 

MTRS0172 

MTRS0173 

MTRS0174 

MTRS0175 

MTRSQ176 

MTRS0178 

MTRSOlgO 

MTRS0181 

MTRS0182 

MTRS0184 

MTRS0185 

MTRS0186 

MTRS0187 

MTRS0188 

MTRS0189 

MTRS0190 

MTRS0191 

MTRS0211 
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3? GOTO (34,36), NX 
34 GOTO 1 6 
36 GOTO ip 
30 GOTO 4(1 
40 GOTO (14,6) , MX 
END 

SUBROUTINE DPMLTX ( A , N A , B , N 8 , C , M , N , K , M A . M 8 , M C ) 
SUBROUTINE 

CALLING SEQUENCE 

CALL DPMLTX ( A,NA»B,N8,C»M,N,K,MA,MB,MC) 


r PBEMLLl ?PL IEP 
t POST MULTIPLIER 
PRODUCT 

NO. ROWS IN A « 
NO. COLUMNS IN A 
NO. COLUMNS IN B 
AND Nfi =. 1 JF A OR 
s ? IF A OR 


MA 8 DIMENSIONED NUMBER of ROWS A 

M8 # dimensioned number of rows b 

MC 8 DIMENSIONED NUMBER OF rows C 

8 NO. ROWS IN c 
A s NO, ROWS IN 8 
B s NO, COLUMNS IN C 
OR 8, RESPECTIVELY, ARE REAL. 

OR 8. RESPECTIVELY, ARE COMPLEX. 


MA,MP,MC ARE SINGLE PRECISION DIMENSIONS, <1/2 O'F ACTUAL CORE 
RFSERV6D FoR REAL DOUBLE PRECISION MATRICES, AND 1/4 OF 
ACTUAL CORF RESERVED FOR COMPLEX DOUBLE PRECfSON MATRICES) 
DIMENSION All), 8(1), CU) 

I AsMC*K 
I BsM A*N*N A 
ICal 

IDsMA*NA 
IHsMC 
I JsMC 
IKsM 
I L»1 
IMsO 

GOTO ( 6, 7 ) , NA 
GOTO (11, 10), NR 
GOTO ( 9 » 6 ) , N B 
IC8? 


GOTO 10 
IHs2*lM 
ICs3 
1 As2*Ia 

DO 30 I s.l , I K , IL 
I NC B I M 

DO 16 J*I,lA,IH 
I Ns I NC 
C( J) = 0 . 

DO 13 Lal.lR.ID 
I N s I N* 1 

C( J)sC( J)+A(L)*B( IN) 
INCsINC+MB 
GOTO (30,18,24) , IC 
IEbI*Ma 
INC«IM 

DO 23 Jsl, I A , I J 
INbINC 
IFbJ+mc 

DO 20 La IE, 18, I D 
I Ns JN+1 
I Gs I N+MB 

C ( I F ) s c ( I F ) ♦ A ( L ) * B ( IN) 
C< J)»C< J)-A(L>*B< 16) 

I NCb I NC + 2*MB 
GOTO 30 
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24 

IE® I *MC 


MTRS0114 


I F s I + m a 


MTRS0115 


I NCs IM 


MTRS0116 


00 ?9 J5 I E > J A , 1 J 


MTRS0117 


I N= INC 


MTRS0H8 

25 

C(J)sO. 


MTRS0120 


DO ?6 UIF, JR, ID 


MTRS0121 


I Ns 1 N*1 


MTRS0122 

2ft 

C(.I)=C( J)*Af L)*B< IN) 


MTRS0123 

29 

INC* I NC*MR 


MTRS0129 

30 

CONTI NUE 


MTRS0131 


RETURN 


MTRS0132 


END 


MTRS0133 


SUrROUTINE SREEPX (HTRUE, u,m, us.fl. 

MODE, N, MD, NC, INDEX, EP) 

MTRS0233 

COMPUTES TRUE MODE AND SWEEPS IT FROM : TWg 

MATRIX. (REAL OR COMPLEX) 

MTRS0223 

HTRUE e TRUE MOCAL COLUMNS# AS COMPUTED. 

U b CYNAMIC MATRIX, 

MTRS0225 


H = SCRIES OF MODIFIED MODAL COLUMNS. 

cl 8 COLUMN OF EIGENVALUES. 

M T R S 0 2 ? 6 


US S SFRIFS of modified modal rows OF 

U. 

MTRS0227 

MODE s MODE NCW rE I NR COMPUTED. N 

s SIZE 

MTRS0228 


MD B DIMENSIONED NUMBER of ROWS of U,US,M,mTRUE 

MTRS0??9 


MX s 1 IF PORBI.eM is real. 


MTRS0230 


8 R IF problem is complex. 


MTRS0231 


DIMENSION H<1), US<1>* J ( It* MTRUE ( 1 ) » 

F l < 1 ) , 6.(4) 

MTHS0235 


msM0DE-1 


MTRS0237 


K1sM*NC*MD 


MTRS 02 38 


no 6 jsI.nc 


MTRS0240 


K*Kl* ( J-l ) *MD 


MJRS0241 


DO 6 L*1 > N 


MTRS0242 


K B K ■# 1 


MTRS0243 

6 

WTRUE(K)BH(K) 


MTRS 0 244 


If ( M ) 31 » 3l » 8 


MTRS0246 

8 

DO 25 Isl.M 


MTRS0247 


U1smc*MD*<M0dE-I ) -NC*MO 


MTRS0248 


GOTO . < 9,11), Nc 


MTRS0249 

9 

GsO, 


MTRS0250 


po 10 Js 1 , N 


MTRS0251 


L*L1+ J 


MTRS0252 


K sK 1 ♦ J 


MTRS0253 

in 

Ge G*US ( L ) **■ TRi.iE ( H ) 


MTRS0254 


GOTO 13 


MTRS0255 

U 

G(l)sO, 


MTRS0256 


G ( 2) bO , 


MTRS0257 


DO 12 J 1 * 1 , N 


MTRS 0 258 


LsLl+jl 


MTRS0259 


KsKl+jl . 


MTRS02*0 


L2SL+MD 


MTRS0261 


K2sK«-MD 


MTRS0262 


G ( 1 ) s G ( 1 ) + L S ( 1. ) *MT RUE < K ) -US < L2 > *HTRUF < K2 ) 

MTRS0263 

1? 

G ( 2 > s G ( 2 ) ♦ L S ( | ) * H T R U fc ( K 2 ) ♦ U S ( L ? ) * H T R U E ( K ) 

MTRS 026 4 

1^ 

KbmODE-I 


MTRS0265 


GOTO ( 14, 19) , NC 


MTRS 0266 

14 

IF (ABSF(FL<K)/fL(M0DE)-1. ) - EP> 15 

i. 15,16 

MTRS0267 

1*5 

.6*1 . 


MTRS0268 


GOTO 17 


MTRS0269 

16 

G=(FL(K)-FL(M0DE) ) / G 


MTRS0270 

1 7 

DO 18 .Is 1, N 


MTRS0271 


K s K 1 * J 


MTRS0272 


L=L1+J 


MTRS0273 

18 

MTRUE ( K ) sh < i ) -G< 1 > • HTRUE ( K ) 


MTRS0274 


GOTO 25 


MTRS0275 

19 

M2*K 


MTRS0276 


Jb2*MODE 


MTRS0277 
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IF ( ABSE( (FL(K*1>*FL(J-1>*FL<K)*FL< J 
1 *EP ) 20,28,22 

2 n IF ( ARSF( (FL(K)*FL( J-l)-FL<K»l)*FLt J 
1 -EP) 21,21,22 

21 G(i)«l. 

G(2)>0, 

GOTO 23 

22 G(3)«G(1)**2*G(2)**2 

G(4)s(FL<K)-FLt J> ) *G<1>« < FL <K-1 > -Ft 
G<1 )*..<< Ft *'K.-1).»FL.< J-l) >*G<1>*<FL(K> 
G<?) = G C 4 ) / 5(3) 

23 DO 24 J1«1»N 

K = K1«-J1 
K2»K* Md 

1 = 1 l + Jl 
L2=L*MD 
G ( 3 ) sHTRUE ( K) 

HTRUE ( K ) * H(L)+ G(2)*HTRUE(K2>- 

HTRUE(K2)= H< L2 ) - G<1>*HTRUE(K2>- 

24 CONTINUE 
2*5 CONTINUE 


) )/(FL(v-l)**2*FL( J)**2>-1. > MTRS0278 

MTRS0279 

>) / < F l. ( J- 1 ) *,* 2 ♦ F L ( J > * * 2 ) ) MTRS028 

MTRS0281 

MTRS0282 

MTRS0283 

MTRS0284 

HTRS0285 

( J- i ) )*G.( 2 ) MTRS0286 

•fl( J) >*G<2> ) / G ( 3 ) MTRS0287 

MTRS0288 

MTRS0290 

MTRS0291 

MTRS0292 

MTRS0293 

MTRS0294 

MTRS0295 

G(1)*RTPUE(K) MTRS0296 

G ( 2) *G ( 3 ) MTRS0294 

MTRS0298 

MTRS0299 


I * 0 

CALL fcPNRMX ( HTRUE ( K 1 + 1 ) , HTRUE (Kl+ 1 ), N, C, I , NO, NC ) 

31 GOTO ( 24 , 32 ) , Nr. 

26 00 29 J* 1 , N 

Ll«( J. 1 )*MD 
1 2 * K 1 + J 
00 29 1 = 1 , N 

L=L 1 +I 

IF ( I « INDEX ) 28 , 27.28 

27 U ( L ) = 0 • 

GOTO 29 

28 K»Kl+ I 
U(L)=U(L)-H(K)*US(L 2 > 

29 CONTINUE 

30 RETURN 

32 00 35 1 = 1 . N 

L 1 =ND*NC*( 1 - 1 ) 

L 2 e K 1 . * I 
J=L 2 *MD 

no 35 Jl=l,N 
L=Ll*Jl 
K 3 =L*N 0 

IF ( Jl- INDEX ) 34 , 33,34 

33 U ( l. > » 0 , 

U ( K 3 ) e 0 • 

GOTO 35 

34 KbkI* J l 
K 2 =K*HD 

U(l )«U<D-H(K)*US(L 2 )*H(K 2 )*US( J) 
U(k 3 )«U(K 3 )«H(K 2 )*US(L 2 )»H(K)*US( J) 

35 CONTINUE 
GOTO 30 

END 

SUBROUTINE MITERS (A, NT APE, N, GUESS, NGUESS, NMODE, VECTOR, 

1 EIGVAL, NITER, NITRSP, EPS P, 

2 US, M.MAXR, NC, A I T KEN , NAKSR,UTRSV ) 

calling sequence 

A s MATRIX, DIMENSIONED (MAXR X 2 *N> - REAL 

(MAXR X 2 *N> - COMPLEX 

N = ORDER OF MATRIX 

GUESSbIST, GUESS VECTOR, DIMENSIONED (MAXR X 1 ) - REAL 


MTRS0301 

MTRS0302 

MTRS0304 

MTRS0305 

MTRS0306 

MTRS0307 

MTRS0308 

MTRS0309 

MTRS0310 

MTRS0311 

MTRS0312 

MTRS0313 

MTRS0314 

MTRS0315 

MTRS0317 

MTRS0319 

MTRS0320 

MTRS0321 

MTRS0322 

MTRS0323 

MTRS0324 

MTRS0325 

MTRS0326 

MTRS0327 

MTRS0328 

MTRS0329 

MTRS0330 

MTRS0331 

MTRS0332 

MTRS0333 

MTRS0334 

MTRS0336 

MTRS0337 

MTRS0375 

MTRS0376 

MTRS0377 

MTRS0341 

MTRS0343 

MTRS0344 

MTRS0346 

MTRS0347 
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c 


c 


c 


(MAXR X 2) - COMPLEX 

NGUESSeO , routine supplies guess vector 

• + 1, GUESS CONTAINS GUESS VECTOR 
NM0DE*NUMBER of eigen SOLUTIONS REQUESTED 
VECTOR* f I GPN VECTORS, DIMENSIONED (MAXR X NMOCE) - REAL 

(MAXR X 2*NM0CE> - COMPLEX 
EIGVAL«F IGtNVALUKS (NMODE X 1) - REAL 

<NMO0F*2 X 1) - COMPLEX 
NITER«NUMBER OF ITERATIONS PER MODE 
NITRSP * MAXIMUM NUMBER Of SINGLE PREC. ITERATIONS 
PPRP a CONVERGENCE CRITERIA FOR SINGLE ROOTS 
US«CMECK EIGENVECTORS, DIMENSIONED (MaXR X NMCDE ) - REAL 

(MAXR X 2*NMCDE> - COMPLEX 

MsViORK I NG AREA OF CORE DIMENSIONED (MAXR X ( A MODE+4 ) - REAL 

(MAXR X 2* ( NM0DE*4 ) - COMPLEX 
WILL CONTAIN CHECK EIGENVALUES, IF REQUESTED 
MAXR » DIMENSIONED NUMBER OF ROWS 
NC * 1, PROBLEM REAL 
s 2, PROBLEM COMPLEX 
AITKFN s AITKEN CONVERGENCE CRITERIA 
NAKSR « NUMBER OF TIMES A I TKEN APPLIED 

DIMENSION A(l). GUESS ( 1 ) , VECTOR(l). EIGVAL(l), NITER(D, USU># 
1 H( 1 ) , NaKSR(1),UTRSV(1) 

define program CONSTANTS and zeros, 

A MODE* 0 

AT*a!TKEN**2 
IF ( EPSP ) 12,9,12 

9 EPSP * ,iE-oe 
1? IF ( NGUESS ) 15,13,15 

13 J1*MAXR*(NC«1) 

DO 14 I»1,N 

K * J 1 ♦ I 
GUESS ( K ) *0 , 

14 GUESS ( I > *1 , , 

15 MODE'* MODE+1 
NAKSR(MODE)*0 
IGOsI 

Nl TER ( MODE >*0 
k1*NC*maXR*( m ODE'»1 ) 
k2*K1+1 

K3*NC*(M0DE-1 )*x 
K4= Nr*MAXR 
K5s k4*NM0DE 
K6*K5*k4 

MOVE FIRST GUESS InTO POSITION 
DO 16 J*l , NC 
JlsMAXR* (U*l ) 

DO 16 1*1, N 

KsKl* J1+ I 
L* Jl* I 

1.6 H(K)sGUESS(L) 

17 N A K * 0 

lfl NITER(MODE)*MTfR(MODE)*1 
nak*nak*1 

I ndex=o 

CALL PPMI.TX (A, NC, H(K2).NCi VECTOR(K2), N,N,1, MAXR, MAXR, MaXR 
CALL NPNRMX (VEcTOR(K2), H(K2) , N, EIGVAL<k 3), INDEX, MAXR, NC 
TEST FOR SINGLE ROOT CONVERGENCE 
DO 23 Jsl.NC 

ji*( j-1)*maxr 

K*K1*J1 

GOTO (24,19,21), NAK 


MTRS0348 
MTRS0379 
MTRS<I350 
MTRS0351 
MTRS0352 
MTRS0353 
MTRS0354 
MTRS0355 
MTRS0356 
MTRS0357 
MTRS0358 
MTRS0360 
MTRS0361 
MTRS0362 
MTRS0363 
MTRS0364 
MTRS0367 
MTRS0368 
MTRS0369 
MTRS0370 
MTRS0371 
MTRS0380 
MTRS0381 
MTRS0408 
MTRSO 4 0 9 
MTRS0411 
MTRS0412 
MTRS0413 
MTRS0417 
MTRS0419 
MTRS0420 
MTRS0421 
MTRS0422 
MTRS0423 
MTRS0425 
MTRS0426 
MTRS0428 
MTRS0429 
MTRS0430 
MTRS0431 
MTRS0432 
MTRS0433 
MTRS0434 
MTRS0435 
MTRS0437 
MTRS0438 
MTRS0439 
MTRS0440 
MTRS0441 
MTRS0442 
MTRS0443 
MTRS0445 
MTRS0446 
MTRS0447 
MTRS0448 
) MTRS0449 
) MTRS0450 
MTRS0452 
MTRS0453 
MTRS0454 
MTRS0455 
MTRS0456 
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19 

L«K5+J1 

MTRS0457 



DO 20 1 * 1 , N 

MTRS0458 



L = l. ♦ 1 

MTRS0459 



K = K*1 

MTRS0460 



IF ( ABSF(M(L)-H(K) ) « EPSP) 20,20,24 

MTR$046l 


2 n 

CONTINUE 

MTRS0462 



GOTO 100 

MTRS0463 


2 l 

DO 2? 1*1, N 

MTRS0464 



KsK*l 

MTRS0465 



IF ( ABSF(US(K)-H(K)> - EPSP ) 22,22,24 

MTRS0466 


2 ? 

CONTINUE 

MTRS0467 


23 

CONTINUE 

MTRS0466 


100 

IF ACCUMULATOR OVERFLOW 108,102 

MTRS0469 


10 ? 

GOTO 56 

MTRS0470 

c 


NO CONVERGENCE, SO TEST MAXIMUM NUMBER OF ITERATIONS. 

MTRS0474 


24 

IF (NtTER(MODE)-NITRSP) 25,46,46 

MTRS0475 

c 


NOT YET EXCEEDED, SO TRY FOR AlTKgNS TIME. 

MTRS0477 


25 

GOTO < 40 , 44 , 31 ) , NAK 

MTRS0478 

c 


TEST FOR AlTKENS CONVERGENCE. 

MTRS0481 


3l. 

GOTO ( 26.36I.NC 

MTRS0482 


26 

DO 2B 1 * 1 , N 

MTRS0484 



J*K5* I 

MTRS0485 



K*K1+l 

MTRS0486 



IF < US(K)*H(J) ) 27,261,27 

MTRS0487 


261 

IF ( HfK)-US(K) ) 32,28,32 

MTRS0488 


27 

IF t { ABSF ( (H(K)-US(K))/ (US(K)-H(J)) )) - A I TKEN ) 

28,26, 32MTRS0489 


2 * 

CONTINUE 

MTRS0490 

c 

all VECTOR ELEMENTS ok, so apply AlTKENS SPfcEEER-UPPER. 

MTRS0492 



DO 3 P I*1,N 

MTRS0493 



J*K5* I 

MTRS0494 



K * K 1 * I 

MTRS0495 



0*(H(K)«2,*US(K)*H( J> ) 

MTRS0496 



IF (G ) 29.30,29 

MTRSO 497 


29 

H ( K ) *H { J) * ( t US < K > - H < J) )**2 / Q) 

MTRS0498 


30 

CONT ! M.JE 

MTRSO 499 



NAKSR(mODE) 8 NAKSR(MODE) ♦ 1 

MTRS0500 



GOTO 17 

MTRS0501 

c 

CONVERGENCE TEST NoT MET. RESTORE AND TRY AGAIN, 

MTRS0503 


3? 

DO 33 L*1»NC 

MTRS0504 



J1*(L"1)*MAXR 

MTRS0505 



DO 33 1*1, N 

MTR$0506 



J*Kl* jl* I 

MTRS0507 



K*K5+J1*I 

MTRS0508 



H(K)mUS< J) 

MTRS0509 


33 

US< J > *M ( J ) 

MTRS0510 



NAK*2 

MTRS0511 



GOTO 18 

MTRS0512 

c 

IF 1 

problem complex, repeat all above fqr complex arithmetic, 

MTRS0514 


36 

DO 38 1*1, N 

MTRS0515 



J*K5 + I 

MTRS0516 



K*K1* I 

MTRS0517 



JJ*J*MAXR 

MTRS0518 



kk*k*maxr 

MTRS0519 



0 * (US(K)-H(J) )**2 ♦ <US(KK)»H( JJ))**2 

MTRS0520 



IF ( 0 ) 37,361,37 

MTRS0521 


361 

IF < (M(K)-US(K) )**2 ♦ (H{KK) - US(KK) )**2 ) 32,38,32 

MTRS0522 


37 

IF ( f (H(K)-nS(K) )**2 ♦ (H(KK) - US ( KK ) ) **2 ) / Q-4T) 38,38,32 MTRS0523 


38 

CONTINUE 

MTRS0525 



DO 39 1*1, N 

MTRS0527 



J=K5+ I 

MTRS0528 



JU*J*MAXR 

MTRS0529 



K = Kl + I 

MTRS0530 
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kk=k*maxR 

a = (H(K>-2,*US(K>*H< J) )**2 * <H<KtO-2,*US<KK)*H< Jj)>**2 
t F ( 0 1 35.39,35 
35 X = H(K) 

N(K)s N(j) - ( ( < (USOO-H( J) )**2»(US(KK)»H( J„ ! ))**2)*(H(K)«2,* 

1 US(K)*H( J) )*(2.*(U5(K)-M< J) )*(US(KK)*M( JJ) )* 

2 < H ( K K ) - 2 . *US < K K ) ♦ H < J J ) ))) / Q) 

M(KK > *H( J.J) • ( ( <?,*(US(K)-H(J))*<US(KK)-H<JJ) ) )*<X-2,* 

1US(K) + H( J) )-( (U«5(K)-W{ J) )**2r <US<KK)-H< JJ)^**2) 

2* ( H ( KK ) *2 , *US ( K< ) *H ( J J > ) ) / 0 ) 

39 CONT I K'llf? 

NAKSR ( MODE ) « NaKSR<HOOE) ♦ 1 
GOTO 17 

40 00 41 Jal , NC 
jlsMAXR* ( J-l ) 

DO 41 I «1 , N 
K*Kl* J1 + I 
L*K5* J1 * I 

41 H(L)«H(K) 

GOTO ( 18.56). 1 GO 

44 DO 45 Jal.NC 
JlsMAXR* ( J-l > 

DO 45 I a 1 » N 
KsKl*Jl*I 

45 US < K ) *R( K > 

GOTO lfl 

C MODE DID NOT CONVERGE IN NORMAL ITERATION 

46 MODE = MODE-1 
PRINT 94, MODE 
GO TO fiO 

56 DO 58 Jsl.NC 

JlsMAXR*! J-lWlNOEX 
DO 58 1*1, N 
K*K1*'MAXR*! J*l>*i 
US ( K ) * A < Jl ) 

5g jiaji*K4 

CALL SWEEPX <V6cT0R,A,H,US»EIGVAL,M0DE,N,MaXR,NC, INDEX.EPSP) 
5g jls(NC-l)*MAXR*TNDEX 
GUESS ( Jl ) * 0 , 

GUESS ( I NDEX ) » 0 , 

6? IF (NMOnE-MODE) 70,70,15 
108 PRINT 131 

MODE * MODE-1 
PRINT 132, MODE 
NMODE a MODE 
GO TO 70 

70 IF ( NTAPE ) 71,80,71 

71 CALL DPMLTX I UTRSV, NC, VECTOR, NO, U5 . N. N. MODE, MAXR, MAXR.MAXR > 

J * 1 

K a 1 

DO 72 I al , MODE 
INDEX e 0 

call npnrmx (ust j>, us <j>.n,h< kj, inoex.maxr.no 

J a J-K4 

72 K a K*NC 

80 PRINT 95 

DO 88 1*1, MODE 

LITRsNITERI I I-NITRSP 
IF { LITR) 81.82,82 

81 LITR*0 
GOTO 85 

82 NITER! TIsNITHSP 


MTRS0531 
MTRS0532 
MTRS0533 
MTRS0534 
MTRS0536 
MTRS0537 
MTRS0538 
MTRS0539 
MTRS0540 
MTRS0541 
MTRS0542 
MTRS0543 
MTRS0544 
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MTRS0549 
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MTRS0551 
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MTRS0553 
MTRS0555 
MTRS0556 
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MTRS0558 
MTRS0559 
MTRS0560 
MTRS0562 
MTRS0564 
M T R S 0 5 6 5 
MTRS0566 
MTRS0615 
MTRS0616 
MTRS0617 
MTRS0618 
MTRS<l6l9 
MTRS0620 
MTRS0621 
MTRS0636 
MTRS0637 
MTRS0638 
MTRS0641 
MTRS0646 
MTRS064 7 
MTRS0648 
MTRS0649 
MTRS0650 


MTRS0666 

MTRS0667 

MTRS0668 

MTRS0669 

MTRS0670 

MTRS0671 

MTRS0672 
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8<5 

03 

84 


GOTO (83,84), Nc 


GOTO 

L»2*l- 


86 

1 


MINT 91 , ( I,6IGVAL( I),MTER( D.NAKSRd)) 


PRINT 96, < I, E I QV At ( L) • E I GV Al.( L+l > » NITER ( I > 


1,NAKSR( I ) ) 

8ft CONTINUE 

IF < NODE > 92,92,88 

88 PRINT 98 
i_sMODE*NC 
CALL MPRINT 
I F ( NT APE ) 

90 PRINT 99 

PRINT 93, ( H < I ) , I * 1 , L ) 

CALL NPRINT (US,N,L,MAXR> 


( VECTOR, n,l,maxr> 

92,92,90 


9? 

RETURN 

93 

FORMAT 

94 

format 

9»5 

1 

format 

96 

format 

97 

format 

9 A 

format 

99 

format 

131 

format 

13? 

format 

END 

END 

end 


(1H ,6Flf,,8) 

( 5H MODE , 1 1 4 , 40H HAS NOT CONVERGED IN MAXIMUM ITERATIONS// 
(1.H17X, MODE 15X, 11H EIGENVALUE 2fX, 

IOHITERaTIONS 11X, 7HAITKENS /> 

<1H 1111, 2619.8, 1122, 1119 > 

(1H 1111, 9X , 1E20.8, 9X , U22, 1119 ) 

( 1 HO / 1H0 46X, 14H EIGENVECTORS ///) 

(1H0 / 1H0 36H CHECK EIGENVALUES AND EIGENVECTORS ) 

(46HlERRnR IN ITERATION SUBROUTINE. . . . ( OVERFLOW )) 

{ 25W CAI CULATION TERMINATED. , H 6, 19H MODES ARE CORRECT.) 
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MTRS0691 

) MTRS0693 
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